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ABSTRACT 

Application  of  classic  triangulation  methods  will  allow 
the  location  of  a  radar  to  be  determined  by  passive  sensors. 
Through  the  use  of  modern  digital  signal  processing  techniques 
this  estimate  can  be  made  in  a  simpler  fashion  using  a 
conventional  receiver. 

In  this  thesis  a  technique  is  developed  for  time 
difference  of  arrival  (TDOA)  estimation  using  a  frequency 
domain  based  correlation  detector  driven  by  an  envelope 
detector.  Time  lag  boundaries  are  defined  on  the  output  of  the 
correlator.  A  fixed  detection  threshold  is  calculated  to 
permit  constant  false  alarm  rate  (CFAR)  detection.  The 
performance  of  the  correlation  detector  is  plotted  as  a 
receiver  operating  characteristic  (ROC)  curve  as  a  function  of 
signal  to  noise  ratio  (SNR) .  An  interactive  MATLAB  software 
program  is  provided  to  perform  either  spectral  domain  or  time 
domain  based  correlation. 

Spectral  domain  based  correlation  uses  the  Fast  Fourier 
Transform  (FFT) .  Implicit  with  the  use  of  the  FFT  are  finite 
arithmetic  internal  processing  errors  which  are  modeled  as 
independent  uncorrelated  noise  sources.  A  method  is  presented 
to  account  for  SNR  degradation  at  the  output  of  the  FFT. 
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I .  INTRODUCTION 

A.  BACKGROUND 

Interception  of  radar  transmissions  by  passive  sensors  is 
one  of  the  responsibilities  of  electronic  intelligence. 
Analysis  of  the  signal  of  a  radar  can  provide  an  estimate  of 
the  distance  and  direction  of  an  emitter  relative  to  an 
intercepting  receiver.  Knowing  the  location  is  an  important 
consideration  in  the  evaluation  of  strategic  and  tactical 
capabilities  of  the  radar. 

Triangulation  is  a  traditional  technique  where  the  angle 
of  arrival  of  a  radar  signal  at  two  spatially  separated 
receivers  provides  an  estimate  of  the  transmitters  location. 
With  the  advent  of  the  computer,  digital  signal  processing 
algorithms  can  be  implemented  that  provide  emitter  range 
information  and  hence  allow  localization  in  a  sequential 
sense.  One  of  these  processing  algorithms  is  called  the  time 
difference  of  arrival  (TDOA)  method. 

B.  TIME  DIFFERENCE  OF  ARRIVAL 

The  TDOA  algorithm  relies  on  two  bistatic  passive  sensors 
attempting  to  receive  the  transmission  of  an  active  radar. 
Reception  by  the  two  sensors  is  preferentially  done  in  the 
main   lobe   of   the   transmitting   radar,   but   due   to   the 


nonsynchronous  nature  of  the  task,  reception  through  one  of 
the  sidelobes  of  the  transmitting  radar  is  expected  with  a 
resulting  decrease  in  received  signal  power. 

Given  a  transmitted  radar  pulse  f(t),  two  physically 
separated  passive  sensors  will  most  likely  receive  this  pulse 
at  times  tx  and  ti+x.  The  time  difference  X  is  proportional  to 
their  differential  distance  relative  to  the  emitter.  This  time 
difference  of  arrival  can  be  used  in  the  estimation  of  the 
differential  range  of  the  emitter  to  the  receivers.  The 
reception  of  the  first  pulse  in  a  pulse  burst  by  the  sensors 
is  paramount  in  performing  the  TDOA  estimate. 

The  concept  of  measuring  the  TDOA  of  a  radar  pulse  between 
two  separated  receivers  can  be  modeled  by  a  hyperbola.  The  two 
sensors  whose  positions  are  known  are  located  at  the  two  focal 
points  of  the  hyperbola.  A  transmitting  radar  is  located  in 
the  area  surrounding  the  two  sensors.  Figure  1  shows  this 
model . 

A  hyperbola  is  the  locus  of  points  in  which  the  distance 
from  any  one  point  on  the  hyperbola  to  one  focus  differs  by  a 
constant  amount  from  the  distance  to  the  other  focus.  This 
differential  distance  is  proportional  to  a  difference  in 
arrival  time  X,  of  a  radar  pulse  as  detected  between  the  two 
sensors.  Once  X  is  measured  the  hyperbola  can  be  drawn.  If  the 
sensors  (moving  aircraft  or  nongeostationary  satellites) 
repeat  the  differential  arrival  time  measurements  at  several 
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time  intervals,  a  series  of  hyperbolas  will  be  generated, 
fixing  the  location  of  the  transmitting  radar.  Alternately, 
the  angle  of  arrival  at  either  sensor  can  be  used  to  define 
the  position  of  the  transmitter  on  the  hyperbola.  A  good 
introduction  into  the  use  of  the  hyperbola  in  navigation  is 
found  in  [Ref.  l:p.  27]. 
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Figure  1.   Time  difference  of  arrival  model. 


This  discussion  assumes  that  the  radar  signal  received  by 
the  two  sensors  is  induced  by  the  direct  wave  from  the 
transmitter.  The  effects  of  refraction  and  reflection  on 
propagating  radar  waves  are  not  covered  here. 


The  purpose  of  this  thesis  is  to  develop  and  test  an 
algorithm  to  estimate  the  differential  arrival  time  x,  of  a 
pulsed  radar  signal  collected  by  two  passive  sensors.  If  the 
pulse  burst  is  collected  in  the  time  domain,  time  correlation 
can  be  used  to  generate  a  TDOA  estimate.  The  Fast  Fourier 
Transform  (FFT)  can  also  be  used  to  perform  the  correlation  of 
the  two  received  signals  (i.e.,  fast  correlation). 

Additive  noise  n(t)  ,  superimposed  onto  a  transmitted  radar 
signal  f (t)  by  the  environment  and  the  radar  receiver,  creates 
a  composite  signal  s (t)  at  the  receiver 

sit)    =  fit)    +  nit)    .  (1) 

The  use  of  the  Fourier  Transform  for  TDOA  estimation  is 
attractive  because  of  the  processing  gain  (PG)  a  radar  signal 
receives  in  being  transformed  from  the  time  domain  to  the 
frequency  domain  during  the  detection  process.  FFT  processing 
gain  for  real  valued  signals  can  be  approximated  by 

PG  idB)    -   (log2  [number  points  in  data  record]  -l)-3  .   (2) 

See  for  example  [Ref.  2:p.  33].  The  signal  to  noise  ratio 
(SNR)  at  the  output  of  a  FFT  (SNRqot.)  is  a  function  of  the  PG 
and  the  SNR  at  the  input  to  the  FFT  (SNRIN)  .  This  relation  is 
given  by 


SNROUT   =  SNRIN   +  PG   .  (3) 

The  corruptive  effects  of  additive  noise  are  reduced  by  the 
transformation,  allowing  the  frequency  domain  detection  of  the 
signal . 

C.   RADAR  TYPES 

The  TDOA  algorithm  is  primarily  applied  to  pulsed  radars. 
The  continuous  wave  (CW)  radar  provides  target  radial  velocity 
through  the  Doppler  shift  of  its  return.  The  CW  radar  does  not 
use  a  pulsed  transmission.  Therefore,  the  TDOA  algorithm 
cannot  not  be  directly  applied  to  the  reception  of  CW  signals. 

Pulsed  radars  can  be  divided  into  two  broad  categories, 
ordinary  low  pulse  repetition  frequency  (PRF)  pulsed  radars 
and  the  pulse  Doppler  radars. 

Low  PRF  radars  can  yield  unambiguous  range  information, 
while  their  radar  returns  do  not  give  target  velocity 
information  directly.  Generally  because  of  the  long  duty  cycle 
of  the  low  PRF  radar,  the  probability  of  receiving  all  radar 
pulses  in  a  pulse  burst  by  two  passive  sensors  is  great.  The 
TDOA  algorithm  can  then  estimate  the  differential  distance  to 
the  emitter  using  the  low  PRF  burst. 

The  pulse  Doppler  radar  uses  coherent  transmission  and 
reception  with  a  moderately  high  PRF  [Ref.  3:p.  17.1].  It 
gives  both  range  and  radial  velocity  information  about  a 
target,  though  not  with  the  same  accuracy  of  the  ordinary 
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pulsed  radar  or  CW  radar  respectively.  In  the  pulse  Doppler 
radar,  as  in  all  non-bistatic  radars,  the  receiver  must  be 
turned  off  while  the  transmitter  is  on.  Because  of  this,  a 
range  ambiguity  exists.  Echoes  having  delay  times  equal  to  an 
integral  multiple  of  the  PRF  will  be  undetected.  Target  echoes 
that  remain  within  this  blanked  time  interval  will  remain 
undetected.  A  second  source  of  range  ambiguity  exists.  For  a 
target  located  at  a  distance  rx  from  the  tracking  radar,  all 
other  targets  in  the  main  beam  of  the  radar  a  distance 

r1(  2rx,   3rx,   4r1#  .  .  . ,  krx 

(4) 

where  k   =  integer     , 

will  appear  to  be  approximately  the  same  distance  from  the 
radar. 

Pulse  Doppler  radars  can  be  divided  into  a  low,  medium  and 
high  PRF  class.  Both  the  low  and  medium  PRF  pulse  Doppler 
radar  transmissions  can  be  passively  intercepted  and  a 
relatively  simple  TDOA  estimate  can  be  formed. 

Reception  of  each  pulse  in  a  pulse  burst  from  a  high  PRF 
pulse  Doppler  radar  by  two  independent  passive  sensors  can  be 
difficult.  If  either  sensor  misses  the  first  pulse,  the  TDOA 
estimate  will  be  in  error.  This  is  true  for  any  radar,  but  has 
a  higher  likelihood  of  occurring  in  the  high  PRF  radars  due  to 
the  shorter  pulse  widths  used.  Should  the  pulse  burst  be  coded 
so  that  a  nonuniform  time  pulse  sequence  is  generated  (i.e., 
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staggered  PRF) ,  then  missed  pulses  by  either  sensor  will  not 
seriously  degrade  the  validity  of  the  TDOA  estimate. 

D.   FAST  FOURIER  TRANSFORM  OUTPUT  SIGNAL  TO  NOISE  RATIO 

The  Fourier  transform  is  a  computationally  intensive 
mathematical  operation.  When  implemented  on  an  FFT  processor 
using  fixed  point  arithmetic,  internal  processing  errors  are 
generated.  These  errors  are  a  function  of  the  transform  size, 
and  the  number  of  bits  used  internally  in  the  FFT  processor. 

The  primary  building  block  of  the  FFT  is  the  butterfly 
algorithm.  There  are  three  dominant  processing  errors  that 
occur  within  the  butterfly.  They  are  scaling,  truncation  and 
trigonometric  errors.  These  can  be  modeled  as  independent, 
uncorrelated  noise  sources.  The  cumulative  effect  of  these 
noise  sources  is  to  reduce  the  potential  SNR  at  the  output  of 
the  butterfly,  and  therefore,  the  FFT  processor. 

This  thesis  will  also  examine  the  effects  of  noise  caused 
by  FFT  processing,  and  will  present  a  method  to  account  for 
the  degradation  of  the  output  SNR  due  to  finite  arithmetic. 
Certain  FFT  implementations  will  keep  the  signal  energy  larger 
with  respect  to  the  three  noise  power  sources  than  others. 


II.  TIME  DOMAIN  BASED  CORRELATION 

A .   BACKGROUND 

The  correlation  function  R(x)  measures  the  degree  of 
similarity  between  two  signals  s(t)  and  g(t)  and  is  described 
by 


1_,T 

2TJ-T 
where  g(x  +   t)    =  git)    displaced  by  the  shift  x 


Rix)    =  limitT  „m  —f*  sit)gix  +  t)  dt 

T         2TJ-T  (5) 


Correlation  can  be  performed  in  a  radar  receiver  between 
a  target  return  corrupted  by  additive  noise,  and  a  replica  of 
the  transmitted  signal  maintained  by  the  radar.  This  replica 
is  designed  into  the  frequency  response  of  the  matched  filter 
of  the  radar  receiver.  [Ref.  4:p.  373]. 

Correlating  two  signals  produces  an  output  that  does  not 
resemble  either  of  the  two  input  signals.  The  shape  of  the 
received  waveform  is  destroyed  during  the  correlation  process. 
The  useful  item  of  correlation  is  an  amplitude  peak,  where  the 
location  along  the  time  axis  of  this  peak  indicates  the  amount 
of  time  that  one  signal  lags  the  other. 

If  the  correlation  process  occurs  between  a  signal  and  a 
delayed  replica  of  itself,  the  process  is  called 
autocorrelation.  If  two  dissimilar  signals  are  correlated,  the 
process  is  called  crosscorrelation. 
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B .   AUTOCORRELATION 

Autocorrelation  is  the  process  by  which  two  identical 
signals,  one  a  delayed  replica  of  the  other,  are  correlated. 
The  autocorrelation  function  of  a  discrete  time  sequence  x(i) 
is  defined  by  [Ref.  5:p.  556] 

W-1-|I| 

^(J)  =  ]T  x(i)x(i   +  I) 

fco  (6) 

where  1  =  shift  operator 

N  =  number  of  data  points  in  x(n) 

This  function  is  not  scaled  with  respect  to  the  shift  operator 
nor  the  number  of  data  points.  Figure  2(a)  shows  the  effect  of 
autocorrelation  on  a  burst  of  13  pulses  of  unit  amplitude. 
Normalized  autocorrelation  produces  a  triangular  envelope 
waveform  of  height  one  at  a  time  lag  of  zero.  The  location  of 
this  peak  indicates  that  the  time  shift  for  maximum  overlap 
between  the  pulse  burst  and  a  shifted  replica  of  itself  is 
zero. 

For  illustration,  a  Barker  coded  sequence  of  pulses  of 
length  13  is  also  autocorrelated  in  Figure  2 (b) .  The  Barker 
code  is  given  by  +++++--++-+-+,  where  +  and  -  could  represent 
0  and  k  radians  carrier  phase  shift,  respectively.  Barker 
coding  of  a  pulse  sequence  constructively  modifies  the 
correlation  output  to  magnify  the  time  lag  peak.  The  zero  lag 
point  has  the  greatest  amplitude.  Adjacent  peaks  are 
significantly  reduced  in  amplitude.  Because  of  this  reduction, 
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Figure  2.  (a)  Normalized  autocorrelation  of  13  equal 
amplitude  pulses,  (b)  Normalized  autocorrelation  of  13  Barker 
coded  pulses. 

if  noise  is  added  to  the  received  signal,  a  Barker  coded  radar 
pulse  burst  will  have  a  correlation  peak  that  is  not  as  easily 
confused  with  adjacent  peaks. 

Now  consider  the  autocorrelation  of  a  random  process.  An 
independent,  identically  distributed  zero  mean  Gaussian 
sequence  with  variance  of  one  is  shown  in  Figure  3.  The 
autocorrelation  of  this  sequence  produces  a  peak  at  zero  time 
lag.  There  is  no  strong  correlation  of  the  Gaussian  noise 
except  at  zero  time  lag.  The  autocorrelation  of  Rayleigh 
distributed  noise  is  shown  in  Figure  4. 
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Figure    3.        (a)    Gaussian   distributed  noise    sequence 
(b)    Normalized   autocorrelation    of   Gaussian   noise. 
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Figure  4.   (a)  Rayleigh  distributed  noise  sequence, 
(b)  Normalized  autocorrelation  of  Rayleigh  noise. 
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Rayleigh  distributed  noise  has  a  mean  (dc)  value  of 


jx 

N  2  (7) 


mean  value  -   o< 

where  o  =  standard  deviation  . 
The  autocorrelation  function  of  Rayleigh  noise  has  two 
components.  A  triangular  envelope  caused  by  the  dc  component 
of  the  Rayleigh  noise,  and  a  weighted  delta  function  at  the 
apex  of  the  triangle  which  indicates  the  time  shift  of  maximum 
overlap  of  the  two  sequences.  Clearly  maximum  overlap  occurs 
at  a  time  shift  of  zero  for  any  autocorrelation  function. 

C .   CROSSCORRELATION 

Crosscorrelation  is  the  process  wherein  two  different 
signals  are  correlated.  The  crosscorrelation  function  of  a 
discrete  time  sequence  is  defined  by 


N-l-UI 

W(I)  =  T    x(i)y(i   +  1) 

i=0 


v._.      ^  (8) 


where  x{i)  and  y{i)  are  two  different  sequences . 
This  function  is  not  scaled  with  respect  to  the  shift  operator 
nor  the  number  of  data  points.  Figure  5  shows  the 
crosscorrelation  of  two  uncorrelated  Gaussian  noise  sequences 
and  two  uncorrelated  Rayleigh  noise  sequences.  Clearly,  no 
dominant  amplitude  peak  is  visible  in  either  crosscorrelation 
plot.  Lack  of  a  dominant  peak  indicates  there  is  no  strong 
correlation  between  the  two  sets  of  noise  sequences. 
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Figure  5.   (a)  Normalized  Gaussian  noise  crosscorrelation. 
(b)  Normalized  Rayieigh  noise  crosscorrelation. 


D.   CORRELATION  COEFFICIENT 

A  measure  of  the  linear  dependence  between  any  two  zero 
mean  stationary  time  functions  x(t)  and  y(t)  is  the 
correlation  coefficient  function  and  is  defined  [Ref .  6:p.  48] 
as 


pxv(x)  = 


£a^2) 


^^m/i^m 


(9) 


where   |p  (t)  |  s  1  for  all   x. 


This  function  is  the  crosscorrelation  function  normalized  by 
the  square  root  of  the  maximum  values  of  the  individual 
autocorrelation  functions.  The  normalizing  factor  is  not  a  lag 
dependent  quantity. 
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Given  two  time  series  x(t)  and  y(t)  where  a  is  any 
positive  number 

pxy(x)  =  lifx(t)  =  a-y(t) 

Pxy(t)  =  -1  if  x(t)    =  -a-y(t)  (10) 

P^(t)  =  0  if  x{t)   and  y(t)   are  uncorrelated. 

A.y 

Another  way  of  defining  a  normalized  correlation 
coefficient  for  a  discrete  time  series  of  finite  duration  is 
given  by 


w-i-UI 


£  x(i)    y(i   +  I) 


i=0 
xy 

'  N-l 


N-l-ll 


(ID 


xE^(i)2\  E  y{i+1) 

>  J=o        >   i=o 


E.   SQUARE-LAW  DETECTION  AND  CORRELATION 

If  x(t)  is  a  real  valued  function  of  time,  the  Hilbert 
transform  of  x(t)  is  defined  by  [Ref.  6:p.  484]  as 


x(t)    =  H  [x(t)] 


J--  7T 


x(u) 


du 


U~U)  (12, 

=  x(t)    *  (-i-) 

71  t 

where  *   =  convolution  operation     . 

The  output  of  a  square-law  envelope  detector  u(t),  can  be 
described  by 
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u(  t)  =  x2(  t)    +  x2(t) 

(13) 

where  x(t)    =  real  valued  time  function 

input  to  square-law  detector     . 

The  correlation  coefficient  function  of  two  time  functions 
x(t)  and  y(t)  is  defined  in  Equation  9.  The  Hilbert  transform 
of  the  correlation  function  is  defined  as 

pxy(T)  =  H     [pxy(T)  ] 

oxoy  (14) 

where  ox  and  oy   =  standard  deviation  of  x  and  y 

respectively     . 

The  correlation  coefficient  for  two  square-law  envelope 
detected  signals,  puv(T),  is  defined  by  [Ref.  6:p.  512]  as 

(15) 

where  u(t)    =  square- law  detected  signal  of  x(t) 
v(t)    =  square- law  detected  signal  of  y(t) 

The  quantity  puv(t)  produces  a  sharper  correlation  peak  than 
pxy(t)  .  This  sharpening  of  the  correlation  function  peak  output 
can  aid  in  locating  the  TDOA  correlation  peak.  An  example  of 
both  pxy(t)  and  puv(T)  is  plotted  in  Figure  6. 

F.   ENVELOPE  CORRELATION  AND  RELATED  STATISTICS 

In  this  thesis,  time  domain  based  correlation  will  be 
performed  on  the  output  of  an  envelope  detector.  An  envelope 
detector  can  be  easily  implemented  in  hardware  using  a  diode. 
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Figure  6.  (a)  Normalized  correlation  coefficient,  (b)  Hilbert 
transform  of  the  normalized  correlation  coefficient  (dotted 
line) . 


To  derive  the  expected  value  and  variance  of  the 
crosscorrelation  function,  two  cases  must  be  considered.  These 
cases  are  noise  only  present,  and  signal  plus  noise  present  in 
both  channels  of  the  correlator.  Under  the  noise  only  case, 
the  statistics  for  the  correlator  output  are  easily  derived. 
We  consider  two  real,  independent  identically  distributed, 
zero  mean  noise  series  at  the  input  to  the  correlator.  Each 
series  has  zero  mean  and  variance  (J2  .  It  is  shown  in  Appendix 
A,  that  the  expected  value  of  the  output  of  the 
crosscorrelation  function  rxy(C)  is 


E[rxy(l)]    =  0 


(16) 
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The  variance  of  the  output  of  the  crosscorrelation  function 
rxv(C)  is 


xy 


02r    =  (7V-|i|)o4 

rxyU)  '   '  (17) 

where  N  =  number  of  data  points 

Since  equal  variances  are  assumed  for  both  time  series,  <J4  is 
the  square  of  the  variance  of  either  noise  series. 

G.   BLOCK  CORRELATION 

Two  time  sequences  can  be  correlated  by  any  of  three 
different  methods.  These  will  be  discussed  below. 

The  first  method  blocks  the  length  of  each  input  sequence. 
Correlation  of  these  blocks  will  provide  an  output  with  a 
variance  given  by  Equation  17.  The  linear  dependence  of  the 
correlation  variance  on  (N- 1  ft  \  )  is  shown  in  Figure  5  for  zero 
mean  Gaussian  noise.  Figure  5  shows  the  correlated  output  at 
C  =0  (zero  time  lag)  is  scaled  by  (N-|0|)=N.  Correlation  output 
at  a  time  lag  of  ±N  are  scaled  by  (N-|n|)=0.  The  linear 
dependance  of  the  correlation  output  variance  on  the  time  lag 
does  not  allow  a  constant  detection  threshold,  which  is  needed 
to  automatically  determine  the  location  of  the  correlation 
peak. 

A  second  method  of  correlation  blocks  a  fixed  length 
sequence  x  (n)  and  correlates  it  with  an  arbitrarily  long 
sequence  y (n) .  The  correlation  output  for  this  algorithm  has 
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a  variance  of  NG4 .  Clearly,  no  scaling  of  the  correlation 
output  variance  as  a  function  of  time  lag  would  occur. 

A  third  method  is  scaling  applied  to  the  output  of  the 
first  technique.  The  correlator  output  is  multiplied  by  a 
weighting  function,  1/ (N- 1  C  |  )  1/2  to  correct  for  the  lag 
dependency  of  the  variance.  The  selection  of  a  constant 
detection  threshold  will  be  addressed  in  Chapter  V. 
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III.  FREQUENCY  DOMAIN  BASED  CORRELATION 

A.   FOURIER  TRANSFORM 

Frequency  domain  techniques  can  be  used  to  perform  fast 
correlation  rather  than  the  time  domain  techniques  previously 
discussed.  The  Fourier  transform  translates  the  discrete  time 
sequence  x (n)  into  the  frequency  domain  . 

For  a  finite  duration  time  series  x (n) ,  the  discrete 
Fourier  transform  (DFT)  is  defined  by  the  pair  of  equations 
[Ref  5:p.  100] 

X(k)    =  £x{n)    e        N       ,    0  ^  k   *  N-l 

n-0 

(18) 

x(n)    =  -  Y  X(k)    e      N       ,0  <.  n  £  N-l   . 
N  &> 

Equation  18  is  used  to  transform  the  N  discrete  sample  values 
into  N  frequency  domain  coefficients. 

The  following  discussion  is  based  on  [Ref.  5:p.  286].  The 
DFT  is  a  computationally  intensive  algorithm  if  implemented 
directly.  The  direct  evaluation  of  the  DFT  requires  N2  complex 
multiplications.  The  time  required  for  the  evaluation  of  a  DFT 
is  proportional  to  N2 . 

Algorithms  have  been  designed  that  take  advantage  of  the 
symmetry  and  periodicity  of  the  DFT,  reducing  the  amount  of 
computation  required  to  solve  the  DFT.  Any  of  these  algorithms 
are  referred  to  as  the  FFT. 
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Generally,  the  FFT  requires  (N/2) log2N  complex 
multiplications.  This  smaller  computational  workload  becomes 
significant  as  N,  the  number  of  data  points  in  the  time 
series,  becomes  larger. 

B.   NORMALIZED  FREQUENCY  DOMAIN  CORRELATION 

In  some  radar  receivers  the  processing  of  pulse  sequences 
occurs  in  the  frequency  domain.  Therefore,  a  frequency 
representation  of  the  time  domain  pulses  may  already  exist. 
The  frequency  coefficients  do  not  necessarily  need  to  be 
transformed  back  into  the  time  domain  to  perform  the 
crosscor relation. 

Frequency  domain  data  can  be  crosscorrelated  by 
multiplying  one  set  of  coefficients  by  the  complex  conjugate 
of  the  second  set  of  coefficients.  The  product  is  translated 
back  to  the  time  domain  through  the  use  of  the  inverse  FFT. 
Frequency  domain  based  normalized  crosscorrelation  is 
described  by 
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w-i       Ji^ni, 


N-l  N-l 

N  £l*wi2  x  £im) 


=  0  otherwise 


N-l 

\Y(k)  |2 

(19) 


where  X(k)  -  FFT  of  the  series  x{n) 

Y{k)  =  FFT  of  the  series  y{n) 

*  =  conjugation 

Pxyil)  -  normalized  correlation  function     . 

To     avoid    circular    correlation,     both    time    domain    pulse 

sequences    are   zero  padded  to  Nz  data  points 

Nz   *  Nx   +  Ny   -  1 

where  iVx  =  number  of  data  points  in  x{n)  '   ' 

iVy  =  number  of  data  points  in  y{n)   . 

To  take  advantage  of  the  processing  speed  of  FFT  algorithms, 
Nz  is  made  a  power  of  two  by  increasing  its  length  with 
additional  zero  padding 

N    =  2m  2.  N    +  N    -   1 
z        *   V   x  (21) 

where  m  is  an  integer  . 

Prior  to  any  zero  padding,  the  dc  component  is  removed  to 
create  zero  mean  pulse  sequences.  Both  sequences,  each  of 
length  Nz  are  then  transformed  to  the  frequency  domain  through 
the  FFT  and  processed  using  the  frequency  domain  correlation 
technique.  A  MATLAB  implementation  of  frequency  domain 
crosscorrelation  is  given  in  Appendix  B  (FreqCorr7 .m) . 
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C.   TIME  DIFFERENCE  OF  ARRIVAL  (FOURIER  DOMAIN) 

The  normalized  crosscorrelation  function  implemented  in 
the  frequency  domain  (Equation  19)  is  used  to  estimate  the 
TDOA.  The  crosscorrelation  of  two  similar  signals  will  produce 
a  dominant  peak,  in  the  output  of  the  correlator.  This  peak  is 
used  as  an  estimate  of  the  TDOA  between  the  two  signals.  For 
two  signals  that  have  not  been  corrupted  by  noise,  the 
location  in  time  of  this  peak  is  easily  determined.  Pulsed 
signals  that  have  been  distorted  by  noise  produce  a  noisy 
correlated  output.  The  output  noise  can  mask  the  TDOA  peak  in 
low  SNR  conditions.  Consequently,  the  probability  of  selecting 
a  noise  sample  instead  of  the  true  TDOA  peak  increases  with 
decreasing  SNR. 

Figures  7,8  and  9  show  the  frequency  domain  based 
crosscorrelation  of  two  pulse  sequences  as  a  function  of 
decreasing  SNR.  A  reception  is  simulated  with  a  pulse  burst 
demodulated  by  one  receiver  fed  into  one  channel,  and  a  pulse 
burst  demodulated  by  a  second  receiver  fed  into  the  second 
channel  of  the  correlator.  Pulse  burst  one  has  a  transmission 
delay  of  zero.  Pulse  burst  two  has  a  delay  of  50  time  units. 
Therefore,  the  total  TDOA  is  50  time  units.  The  TDOA  peak  is 
clearly  seen  at  the  high  SNR  of  +7dB  and  +3dB.  At  OdB  SNR, 
noise  peaks  exceed  in  amplitude  the  true  TDOA  peak  causing  an 
incorrect  estimate  if  one  was  selected. 
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Figure  7.  Normalized  frequency  domain  based  crosscorrelation 
SNR  +7dB. 
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Figure  8.  Normalized  frequency  domain  based  crosscorrelation. 
SNR  +3dB. 


Visually  estimating  the  location  of  the  peak  with  the 
maximum  amplitude  in  the  output  of  the  correlator  is  a  easy 
task.  The  decision  rule  states  that  the  point  with  the  maximum 
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Figure  9.  Normalized  frequency  domain  based  crosscorrelation 
SNR  OdB. 


amplitude  is  selected  as  the  TDOA  estimate,  provided  the 
amplitude  exceeds  a  predetermined  threshold. 

When  many  blocks  of  data  must  be  processed,  visual 
inspection  of  the  correlation  output  is  not  possible.  Chapter 
V  discusses  the  design  of  a  threshold  for  automated  data 
processing. 
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IV.  NOISE  AND  SIGNALS  IN  RECEIVER 

A.   NOISE  AND  SIGNALS  IN  THE  RADAR  DETECTOR 

The  purpose  of  a  radar  receiver  is  to  process  and  detect 
the  reflected  energy  from  a  target  being  tracked  by  the  radar. 
An  I/Q  demodulator  is  assumed  in  the  radar  receiver,  located 
after  the  last  intermediate  frequency  (IF)  amplifier  and  prior 
to  the  video  signal  display.  Data  to  be  processed  for  TDOA  can 
be  obtained  from  three  locations  on  the  I/Q  demodulator.  These 
locations  include  the  envelope  output,  the  envelope  squared 
output  and  the  I/Q  channel  data. 

White  Gaussian  distributed  noise  is  assumed  to  be  added  to 
radar  pulses  during  their  transmission  through  the  atmosphere 
and  reception  by  the  receiver.  This  noise  contributes  to  the 
degradation  of  the  SNR  and  is  caused  by  fluctuations  in  the 
target  radar  cross-section  and  the  loss  of  signal  energy  due 
to  the  propagation  distance  between  the  target  and  the 
receiver. 

At  typical  radar  frequencies  the  front  end  amplifier  of 
the  radar  receiver  contributes  to  this  noise  power.  The  noise 
power  for  a  thermally  noise  limited  receiver  is  defined  as  N0, 
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N0  =  k  Tsys  B 

where  k  =  Boltzman's  constant 
B  -  IF  bandwidth 

»        -  (22) 

T        =  T      +T 

rar  =  receive  antenna  temperature 
Te   =  effective  noise  temperature 
of  the  receiver     . 

The  total  noise  is  modeled  as  zero  mean  Gaussian  distributed 
noise  with  a  probability  density  function  (pdf)  of 

1    ^  (23) 

\/2-7T-o2 

where  a2  is  the  noise  variance  due  to  all  noise  sources  (i.e., 
receiver  front  end,  transmission  noise) . 

For  an  operational  radar,  either  a  pulsed  signal  plus 
Gaussian  noise  or  Gaussian  distributed  noise  alone  is  assumed 
at  the  input  of  the  I/Q  demodulator.  The  probability 
distribution  of  the  demodulated  signal  is  dependant  on  which 
of  the  three  outputs  of  the  I/Q  demodulator  is  used. 

1 .   Coherent  detection 

If  the  exact  frequency  and  phase  of  the  received  pulse 
burst  is  known,  and  matches  the  frequency  and  phase  of  the 
local  oscillator  in  the  receiver,  then  the  signal  can  be 
coherently  detected.  Coherent  detection  simplifies  the 
probability  description  of  the  signal  and  noise  at  the  output 
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of  the  I  and  Q  channels.  The  I  and  Q  sequences  have  an 
independent  Gaussian  distribution,  and  can  be  described  as 

1,0   ~  N(  {m^n)  ,mQ{n)  }  ,o2)  (24) 

where  nij^in)  ,   mQ{n)    =  I  and  Q  mean  values  respectively 

o2  =  noise  variance  due  to  all  noise  sources  . 

Both  the  I  or  Q  channel  can  be  processed  using  the  correlation 
algorithm  followed  by  threshold  detection. 

Due  to  the  passive  nature  of  signal  reception  and 
demodulation,  the  exact  frequency  and  phase  of  the  received 
pulse  burst  are  not  available  at  the  receiver.  Therefore, 
coherent  detection  is  not  feasible  for  TDOA  estimation. 
2 .   Envelope  Squared  Detection 

For  a  zero  mean  Gaussian  distributed  noise  input  to 
the  I/Q  demodulator,  I2  and  Q2  each  have  a  chi-squared 
distribution  with  one  degree  of  freedom.  Their  pdf ' s  can  be 
described  by 

PY(y)    = — e-y/2°2  foryzO 

c/2ny 

-   0  for  y   <  0        (25) 

where  o2  =  variance  of  I  and  Q 
y   =  I2  or  Q2     . 

The  envelope  squared  output  is  the  sum  of  I2  and  Q2  and  is 

therefore  chi-squared  with  two  degrees  of  freedom.  This  pdf 

can  be  described  by 
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e-z/202 

pz{z)    =   for  z  -2.   0 

2°2  (26) 

where  z   =  J2  +  Q2     . 

Envelope  squared  detection  is  attractive  because  the  Hilbert 
transform  can  be  used  to  enhance  the  TDOA  peak  in  the 
correlator's  output.  Unfortunately  in  a  low  SNR  environment 
the  noise  power  will  also  increase  as  the  noise  contribution 
is  squared. 

3 .   Envelope  Detection 

A  common  form  of  receiver  demodulation  at  IF 
frequencies  is  the  envelope  detector.  For  all  simulations  in 
this  thesis,  envelope  detection  will  be  assumed.  The  output  of 
the  envelope  detector  is  the  modulation  envelope  of  the 
received  pulse  burst.  For  Gaussian  distributed  noise  at  the 
input  to  the  I/Q  demodulator,  the  noise  at  the  envelope  output 
is  Rayleigh  distributed  with  pdf  given  by 


i2 


pa{u)    =  -^e  20  ,  u   >  0 

-  0,        otherwise  (27) 


where  a2   =  variance  of  Gaussian  noise 

u   =  sfz   =  s/l2   +  Q2     . 

For  a  pulsed  sinusoid  plus  Gaussian  noise  at  the  input  to  the 
I/Q  demodulator,  the  pdf  of  the  envelope  at  the  output  is 
Rician  distributed.  The  Rician  pdf  is  given  by 
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lA 
a2'  "°o2 


pR(r)    =  4e      2°2      I0(^)    ,    r  2  0 
o2  o2 

=  0,  otherwise  (28) 

where  A  =  amplitude  of  received  sinusoid 
o2  =  variance  of  Gaussian  noise 


where  lAz)    a  —    P*  ez  cos(  e  >  d6 
0      2n  Jo 


(29) 


and  I0(z)    =  modified  Bessel  function  of  the 
first  kind  of  zero  order     . 


Figure  10  shows  a  block  diagram  of  the  I/Q  demodulator  and  the 
various  pdfs  in  the  radar  receiver. 

B.   ENVELOPE  CORRELATION  OF  RAYLEIGH  NOISE 

Rayleigh  distributed  noise  is  produced  at  the  output  of 
the  envelope  detector  under  noise  only  conditions.  This  noise 
is  then  processed  by  the  correlator  during  the  TDOA  estimate. 
In  all  follow  on  discussions  and  simulations,  the  envelope 
processing  scheme  is  used. 

The  autocorrelation  of  Rayleigh  noise  as  discussed  earlier 
produces  a  large  triangular  bias  in  the  output  which  decreases 
as  the  time  lag  increases.  This  bias  is  produced  by  the  dc 
component  of  the  Rayleigh  distributed  noise. 

This  triangular  function  contains  no  information  relative 
to  the  TDOA  estimate.  The  large  numerical  values  in  the  bias 
can  potentially  limit  the  dynamic  range  in  the  FFT  based 
correlation.   A   constant,   detection   threshold   cannot   be 
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Figure  10.   I/Q  demodulator  in  the  radar  receiver. 


implemented  at  the  output  of  the  correlator  with  this 
triangular  function  present.  For  these  three  reasons,  it  is 
advantageous  to  force  the  Rayleigh  or  Rician  distribution  to 
have  a  zero  mean.  This  will  remove  the  triangular  bias  formed 
by  the  correlation  process.  This  modification  of  the  received 
pulse  burst  does  not  degrade  the  TDOA  estimate.  For  all 
realizations  (i.e.,  blocks  of  data),  the  dc  component  of  the 
Rayleigh  noise  will  be  removed. 
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1 .   Shifted  Rayleigh  Noise 

Forming  a  Rayleigh  distributed  time  series  x  (n)  and 
subtracting  its  expected  value  produces  a  shifted  Rayleigh 
sequence.  This  algorithm  is  described  by 


x1  (n)    =  x{n)    -  o* 
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(30) 


where  x1(n)   =  shifted  Rayleigh  noise  sequence 

o  =  Gaussian  noise  standard  deviation 


The  shifted  Rayleigh  time  series  has  a  zero  mean.  The  new  time 
series  and  the  corresponding  time  domain  based  autocorrelation 
function  are  shown  in  Figure  11. 
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Figure  11.  (a)  Shifted  Rayleigh  time  series,  (b)  Normalized 
shifted  Rayleigh  time  domain  based  autocorrelation. 


Histograms  of  a  Rayleigh  pdf  and  a  shifted  Rayleigh 
pdf  are  shown  in  Figure  12.  Clearly,  the  effect  of  removing 
the  dc  term  from  the  Rayleigh  distributed  time  series  creates 
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(a)  Rayleigh  pdf.  (b)  Shifted  Rayleigh  pdf. 


an  autocorrelated  output  that  does  not  have  the  limiting 

triangular  bias.  The  impulse  at  zero  time  lag  dominates  the 

correlation  output.  A  constant  detection  threshold,  at  least 

over  a  small  number  of  range  bins,  could  be  implemented  with 

this  type  of  correlation  output. 

2 .   Frequency  Domain  Based  Correlation  of  Shifted  Rayleigh 
Noise 

The  frequency  domain  based  autocorrelation  of  xx(n)  is 

now  examined.  Sequence  xx  (n)  is  zero  padded  as  described  in 

Equation  21.  A  frequency  domain  series  is  created  using  the 

FFT  as  described  by  Equation  18.  The  autocorrelation  estimate 

is  derived  using  Equation  19  and  is  plotted  in  Figure  13.  The 

results  of  frequency  domain  based  autocorrelation  are  very 

similar  to  the  results  of  time  domain  based  autocorrelation 

using  shifted  Rayleigh  noise. 
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Figure  13.  Normalized  frequency  domain  based  autocorrelation. 

3 .   Spectral  Interpolation 

With  the  incorporation  of  the  FFT  in  modern  radar 
receivers,  it  is  possible  that  the  received  pulse  burst  has 
been  transformed  into  frequency  data  as  a  byproduct  of 
detection.  Processing  time  would  be  lost  in  converting  this 
data  back  into  a  time  series  to  be  correlated.  Of  course, 
there  is  no  guarantee  that  the  dc  component  has  been  removed. 
In  this  section  an  algorithm  is  given  to  perform  frequency 
domain  based  correlation  on  a  Rayleigh  distributed  noise 
sequence  while  avoiding  the  triangular  correlation  bias 
created  by  the  (possible)  dc  component. 

A  Rayleigh  distributed  noise  series  x(n),  N  points 
long,  is  transformed  into  an  ordered  sequence  of  N  complex 
coefficients  using  the  FFT  algorithm.  In  general,  a  sequence 
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X(k)  of  coefficients  is  indexed  from  zero  to  N-l  and  is 
described  by 

X(k)    =  FFT  {   x(n)    }  ,      0  £  k  zN-l 

(31) 

where  k  =  frequency  index     . 

The  FFT  creates  the  zero  indexed  term  X(0)  by  summing  all  the 
elements  in  x (n) .  X(0)  is  real  and  can  be  viewed  as  the  dc 
component  of  x (n) .  By  removing  the  X(0)  coefficient  from  the 
frequency  sequence  we  have  in  effect  removed  the  dc  bias  from 
the  time  series  x  (n)  .  At  this  point,  the  modified  Fourier 
transform  X(k)  approximates  the  DFT  of  a  shifted  Rayleigh 
series . 

Time  domain  correlation  of  an  N  point  series  produces 
2N  data  points.  Frequency  domain  based  correlation  of  an  N 
point  sequence  produces  N  data  points.  To  force  frequency 
domain  based  correlation  to  equal  time  domain  based 
correlation  (i.e.,  avoid  circular  correlation),  the  number  of 
terms  in  the  frequency  series  must  be  doubled.  The  frequency 
series  is  expanded  by  storing  each  pair  of  complex  frequency 
coefficients  at  a  location  twice  the  value  of  the  original 
index. 

As  a  first  approximation,  the  data  points  between 
adjacent  coefficients  are  linearly  interpolated  in  the 
expanded  array.   Figure  14  illustrates  the  effect  of  FFT 


34 


interpolation  and  the  zeroing  of  the  dc  terms.  Correlation 
results  as  described  by  Equation  19  are  shown  in  Figure  15. 
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Figure  14.     FFT   interpolation,   and   zeroing  of  the   dc 
coefficient. 
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Clearly,  the  triangular  bias  in  the  correlation  output  has 
been  removed.  A  small  impulse  at  a  time  lag  of  -N  is  observed, 
which  can  be  disregarded. 

For  block  processing,  typically  only  correlation 
products  within  ±10%  of  the  zero  time  lag  position  are 
considered  to  be  statistically  reliable.  Therefore,  the  TDOA 
estimate  is  formed  within  this  correlation  band.  The  impulse 
at  -N  falls  outside  of  the  TDOA  evaluation  area  and  will  not 
cause  a  false  detection. 

C.   CORRELATION  OF  PULSES  IN  ADDITIVE  RAYLEIGH  NOISE 

Time  domain  based  correlation  of  a  pulse  train  in  additive 
Rayleigh  noise  should  produce  the  same  result  as  frequency 
domain  based  interpolation  followed  by  an  equivalent 
correlation.  A  comparison  was  performed  using  a  50%  duty  cycle 
pulse  train  in  additive  shifted  Rayleigh  noise. 

Time  correlation  was  performed  on  the  pulse  burst.  The  dc 
component  of  the  noisy  pulse  train  is  removed  forming  a  zero 
mean  signal.  The  time  domain  based  autocorrelation  of  this 
sequence  is  shown  in  Figure  16(a) . 

The  pulse  burst  was  also  transformed  into  the  frequency 
domain  using  the  FFT  with  the  dc  component  intact.  N  Fourier 
coefficients  are  expanded  and  the  X(0),  X(±l)  coefficients 
zeroed  out.  The  sequence  is  correlated  after  interpolation 
using  Equation  19  and  plotted  in  Figure  16(b) . 
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Figure  16.  (a)  Normalized  time  domain  based  autocorrelation, 
(b)  Normalized  frequency  domain  based  interpolation  and 
autocorrelation.  Zero  time  lag  between  each  channel. 

Clearly,  the  two  correlation  algorithms  produce  similar 
outputs.  Both  have  a  dominant  peak  in  their  output  for  TDOA 
estimation.  Frequency  domain  based  correlation  displays  a 
exponential  decay  of  the  sidelobes  about  the  main  peak.  This 
effect  is  not  observed  in  time  autocorrelation,  and  is  caused 
by  the  first  order  interpolation  used  in  the  frequency  domain 
based  method.  A  more  robust  interpolation  algorithm  will 
eliminate  the  exponential  decay,  yielding  the  desired  linear 
decay. 

Next,  time  and  frequency  domain  based  crosscorrelation  are 
performed  on  two  pulse  sequences.  One  pulse  sequence  lags  the 
other  by  20  units.  Both  time  and  frequency  domain  based 
crosscorrelation  of  the  two  sequences  are  plotted  in  Figure 
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17.  Clearly,  both  correlation  algorithms  produce  a  peak  at  20 
units.  The  nonlinear  decay  of  the  sidelobes  surrounding  the 
dominant  peak  is  again  seen  in  frequency  domain  correlation. 
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Figure  17.  (a)  Normalized  time  domain  based  crosscorrelation. 
(b)  Normalized  frequency  domain  based  interpolation  and 
crosscorrelation.  20  unit  time  lag  between  each  channel. 

In  summary,  the  frequency  domain  based  interpolation  and 
crosscorrelation  algorithm  produces  an  output  comparable  to 
time  domain  based  crosscorrelation.  It  offers  a  method  to 
correlate  pulse  bursts  collected  in  the  frequency  domain. 
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V.  THRESHOLD  DESIGN  FOR  ENVELOPE  CORRELATION 

A.  INTRODUCTION 

A  threshold  at  the  output  of  a  correlator  is  used  to 
estimate  the  location  of  a  TDOA  peak  in  a  constant  false  alarm 
rate  (CFAR)  receiver.  The  value  of  the  threshold  is  a  function 
of  the  correlation  output  mean  and  variance. 

B.  CORRELATION  OUTPUT  AND  CONSTANT  VARIANCE 

If  the  variance  of  the  correlator,  when  driven  by  two 
signals  at  the  input,  maintains  a  fixed  value  about  some  mean, 
a  constant  threshold  above  the  mean  can  be  used.  This  allows 
the  design  of  a  constant  false  alarm  rate  detector.  The 
largest  correlation  peak  that  also  crosses  the  threshold  is 
defined  as  the  TDOA  estimate.  The  constant  threshold  design  is 
the  easiest  to  implement. 

C.  CORRELATION  OUTPUT  AND  TRIANGULARLY  SHAPED  VARIANCE 

The  crosscorrelation  of  two  zero  mean  noise  sequences  with 
the  same  number  of  data  points,  produces  an  output  whose 
variance  has  triangular  form.  Hence,  a  constant  detection 
threshold  is  not  easily  implemented.  Instead,  a  threshold  must 
be  designed  that  has  a  slope  that  matches  the  slope  of  the 
correlation  output.   The  design  of  a  sloping  threshold  is 
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difficult,  particularly  if  the  slope  of  the  correlation  peak 
changes  with  time  (i.e.,  changing  noise  environment). 

A  solution  is  to  mathematically  force  the  triangular 
correlation  variance  to  be  flat.  This  is  done  by  using  a 
bowtie  correction  of  the  form  1/ (N- | C  |  )  1/2 .  A  bowtie  correction 
normalizes  the  correlation  output  to  have  a  constant  (lag 
independent)  variance.  This  design  suffers  from  time  changes 
in  the  correlation  output  forcing  the  bowtie  correction  to  be 
recalculated.  Note  that  the  amplitude  of  the  TDOA  peak  will 
itself  be  dependant  on  the  type  of  correction  to  the 
correlation  output. 

D.   THRESHOLD  BASED  ON  ZERO  LAG 

The  zero  lag  threshold  method  uses  the  variance  of  the 
zero  lag  correlation  output  to  calculate  a  constant  detection 
threshold.  This  is  indicated  in  Figure  18. 

The  correlation  output  becomes  statistically  less  reliable 
the  further  removed  the  correlation  products  are  from  the  zero 
time  lag  position.  To  use  correlation  intelligently,  the 
correlation  output  is  typically  examined  ±10%  of  the  distance 
from  the  zero  time  lag  position.  Consequently,  only 
correlation  output  ±0 . IN  from  zero  time  lag  will  be  used  in 
TDOA  estimate  simulations. 

The  detection  method  in  this  thesis  uses  a  constant 
threshold  against  the  correlation  output.  The  fixed  threshold 
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Figure  18.   Correlation  output  with  zero  lag  threshold. 

was  selected  because  the  correlation  output  variance  is,  to  a 
first  approximation,  nominally  flat  within  the  ±0.1N  boundary. 
The  detection  scenario  that  drives  this  thesis  defines  a 
minimum  radar  pulse  width  of  one  microsecond  (p.sec)  and  a 
burst  length  of  four  millisecond  (msec)  .  A  50%  duty  cycle 
radar  pulse  is  assumed.  There  are  2000  pulses  per  burst,  which 
when  correlated,  will  generate  an  output  that  is  8000|isec 
(±4000|isec)  in  length.  If  ±10%  of  the  correlation  output 
around  zero  time  lag  are  compared  to  the  threshold  value, 
±400|xsec  will  be  tested.  For  radar  pulses  traveling  at  the 
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speed  of  light,  a  maximum  distance  of  ±120km  is  therefore 
examined  to  form  a  TDOA  estimate. 

Correlation  values  at  lags  -0.1N  to  0 . IN  will  be  compared 
in  magnitude  to  the  zero  time  lag  based  calculated  threshold 
value. 

E.   THRESHOLD  CALCULATION  FOR  CFAR 

The  detection  threshold  is  calculated  using  the  output 
terms  of  the  crosscorrelation  algorithm  when  noise  only  is 
injected  into  the  correlator.  For  convenience,  Equation  8  is 
repeated  here, 


AT-m-i 

E 

i  =  0 


rxy(l)  ]T  x±yi+1 

(32) 


x±  and  yi+1    are  zero  mean  Rayleigh   {shifted) 
distributed  noise  sequences     . 

The  correlation  function  at  all  lags  of  interest  is  thought  to 
be  Gaussian  distributed.  According  to  Equation  32,  every 
estimate  is  the  sum  of  a  fairly  large  number  of  products 
(i.e.,  central  limit  theorem) .  The  average  value  of  the 
correlation  output  between  the  ±0.1N  endpoints  is 


(33) 


where  a   =  number  of  data  points  between  the 
0  and  0  .  IN  lag  posi  tions  in  the 
correlation  output     . 


The  variance  of  rxy(l)  is  calculated  from 

42 


°2t  =  — - —  y  ^xvd 

2a  +   l  £?a       xy 


)  -  t  )2   .  (34) 


For  Gaussian  distributed  noise  the  probability  of  false  alarm 
(Pfo)  is  a  function  of  the  detection  threshold 

Pfa   =    1      fe^  dx 

v/2¥otJ^  (35) 

where  T  -  detection  threshold    . 
For  a  given  Pfa/.  the  threshold  can  be  calculated  from 

T  =  erfc'1   (Pfa)at 
where  erfc  =  complementary  error  function     . 


(36) 


Equation  36  relates  the  threshold  to  the  Gaussian  distributed 
noise  variance. 

F.   VERIFICATION  OF  CFAR  THRESHOLD  AND  ROC  CURVE 

Using  the  above  equations,  thresholds  were  calculated  for 
the  time  correlation  detector  given  in  Equation  32  for  three 
different  Pfa's.  A  MATLAB  simulation  was  performed  to  compare 
the  measured  false  alarm  rate  of  the  time  correlation  detector 
to  the  theoretical  false  alarm  rate. 

1 .   Pf,  Simulation 

Two  100  point  sequences  of  zero  mean  Rayleigh  noise 
were  injected  into  the  correlation  detector.  The  output  of  the 
correlation  detector  between  the  -0.1N  and  +0.1N  lag  points 
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was  compared  to  the  calculated  threshold  to  measure  the  false 
alarm  rate.  The  number  of  points  in  the  correlation  detector 
output  that  exceeded  the  threshold  for  a  given  Pfa  were 
recorded.  The  experimental  P£a  is  given  as  the  ratio  of 
improper  threshold  crossings  to  the  total  number  of  monitored 
range  cells  of  the  correlator  output  (see  Table  1) .  For  each 
Pfa,  500  realizations  were  performed.  Both  theoretical  and 
experimental  results  are  listed. 


TABLE  1 
THEORETICAL  PROBABILITY  OF  FALSE  ALARM 

Simulation  # 

Pf.=  0.01 

Pfa=0.001 

Pfa=0.0001 

1 

0.0100 

0.00076 

0 

2 

0.0099 

0.00090 

0 

3 

0.0110 

0.00085 

0 

The  data  in  Table  1  show  that  a  detection  threshold 
can  be  established  for  a  correlation  detector,  using  the 
assumption  of  a  Gaussian  distributed  output  for  Pfa'  s  of  0.01 
and  0.001.  For  the  lower  Pfa  value  of  0.0001,  the  threshold 
calculation  fails,  but  in  a  positive  sense.  The  measured  Pfa 
of  0  is  below  the  theoretical  Pfa  of  0.0001.  Either  the 
correlation  variance  is  smaller  than  Equation  36  predicts,  or 
more  terms  must  be  summed  in  the  correlation  output  to  better 
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approximate  the  Gaussian  pdf  as  indicated  by  the  central  limit 
theorem. 

2 .   TDOA  Simulation 

The  output  of  the  envelope  detector  will  either  be 
Rayleigh  or  Rician  (non-central  Rayleigh)  distributed.  In 
those  portions  of  the  data  where  there  is  no  signal  energy, 
the  output  will  be  Rayleigh  distributed.  Where  the  signal  is 
present,  the  envelope  detector  output  will  be  Rician 
distributed. 

The  purpose  of  the  time  domain  based  correlation 
detector  is  to  measure  the  TDOA  of  pulsed  sequences  embedded 
in  Rayleigh  noise.  Noise  in  low  SNR  environments  will  mask  the 
TDOA  peak.  A  MATLAB  simulation  was  performed  to  examine  the 
performance  of  the  correlation  detector  in  a  decreasing  SNR 
environment . 

Two  pulse  bursts  embedded  in  additive  Rayleigh  noise 
were  created.  The  second  pulse  burst  lagged  the  first  by  three 
time  units.  Each  pulse  burst  is  a  sequence  of  ten  pulses,  each 
having  a  period  of  ten  units.  Each  pulse  in  the  burst  had  a 
50%  duty  cycle. 

The  SNR  established  the  amplitude  relationship  between 
the  noise  and  signal  sequences.  First,  two  sequences  of  zero 
mean  unit  variance  Gaussian  distributed  noise  were  created. 
Rayleigh  noise  was  generated  from  the  two  Gaussian  sequences. 
The  Rayleigh  noise  had  a  mean  of  1.2533  and  a  variance  of 
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0.4292.  The  theoretical  second  moment  of  the  Rayleigh  noise  is 
two  and  was  verified  in  the  MATLAB  simulation.  The  second 
moment  is  defined  as  the  average  power  of  the  Rayleigh  noise. 
The  amplitude  of  the  nonzero  elements  in  each  pulse  burst  was 
then  scaled  to  meet  the  desired  SNR,  using  the  measured 
average  power  of  the  Rayleigh  noise. 

A  detection  threshold  was  calculated  based  on  shifted 
Rayleigh  noise  using  Equation  36  where  a2  was  assumed  to  be 
known  (i.e.,  simulation  information).  In  practical 
applications,  a2  would  be  estimated  from  the  data.  The  SNR  of 
the  two  bursts  were  varied  from  OdB  to  18dB  in  the  different 
simulations.  For  each  SNR,  the  two  sequences  were  correlated 
using  Equation  32. 

Graphs  of  the  correlation  output  for  a  SNR  of  ldB  are 
given  in  Appendix  C  for  Pfa'  s  of  0.01,  0.001,  0.0001  and 
0.00001.  The  threshold,  which  is  a  function  of  the  Pfa  and 
noise  variance,  is  shown  to  increase  with  decreasing  Pfa. 

The  correlation  peak  at  a  time  lag  of  three  was  the 
correct  TDOA  estimate.  If  the  maximum  correlation  peak  in  the 
output  was  the  third  lag  point,  a  correct  TDOA  estimate  was 
obtained.  For  each  SNR,  30  TDOA  estimates  were  made.  The 
correlation  function  output  100  data  (lag)  points  per  TDOA 
estimate.  The  TDOA  estimate  was  made  on  ±10  data  points  on 
each  side  of  the  zero  time  lag  position  in  the  correlation 
output.  The  percent  correct  TDOA  estimate  for  each  SNR  was 
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calculated  and   is   presented   as   a   receiver   operating 
characteristic  (ROC)  curve  in  Figure  19. 

The  ROC  curve  shows  that  as  the  SNR  increases,  the 

performance  of   the   correlation   detector   increases.   The 

performance  of  the  correlation  detector  in  Figure  19  can  be 

improved  if  multiple  sequential  looks  at  the  same  SNR  were 
allowed. 
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Figure  19.   Percentage  of  estimated  TDOA  correct  versus 
SNR  for  time  domain  based  correlation  detection. 
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VI.  FFT  ERROR  ANALYSIS 

A.  INTRODUCTION 

This  chapter  discusses  several  sources  of  processing 
errors  occurring  within  the  FFT,  and  the  effect  of  these 
errors  on  the  FFT  output  signal-to-noise  ratio. 

The  computation  of  a  FFT  using  fixed  point  arithmetic, 
generates  internal  processing  errors.  These  errors  are  due  to 
the  fixed  number  of  bits  in  the  trigonometric  look  up  table, 
the  truncation  of  arithmetic  results  during  addition  and  the 
scaling  of  results  during  the  individual  butterfly 
calculations . 

The  error  sources  can  be  modeled  as  independent  additive 
noise  components  at  the  output  of  the  FFT.  Their  net  effect  is 
to  reduce  the  ideal  SNR  at  the  output  of  the  FFT. 

Each  of  these  processing  error  sources  will  be  described 
below.  The  material  in  this  chapter  borrows  from  [Ref .  7] . 

B.  COSINE/SINE  TABLE  NOISE 

The  computation  of  a  butterfly  algorithm  in  the  FFT 
involves  multiplication  by  the  complex  coefficient  (twiddle 
factor) 
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A  FFT  butterfly  is  shown  in  Figure  20. 
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Figure  20.   FFT  butterfly  calculation. 

These  trigonometric  coefficients  are  precalculated,  quantized 
to  (B)  bits  (where  B  is  the  size  of  the  computer  word 
including  sign  bit) ,  and  stored  in  memory  for  future  table 
lookup  during  the  FFT  processing.  The  coefficients  W°  =  1  and 
W°/4  =  -j  can  be  obtained  in  an  exact  form  and  therefore  do  not 
contribute  to  the  quantization  noise.  Quantization  of  the 
other  (twiddle)  coefficients  stored  in  the  sine  and  cosine 
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table  introduce  noise  into  the  output  of  the  FFT.  The  noise 
power  created  by  truncation  of  the  sine  and  cosine 
coefficients  is  defined  by  [Ref.  7]  to  be 


a2w   =  ±2-2B 
3 

where  o2w  =  trigonometric  noise  power 

B  =  number  bits  in  computer  word 


(38) 


C.   TRUNCATION  NOISE 

The  DFT  of  a  finite  duration  sequence  {x(n)  }  is  defined  in 
Equation  18  and  is  repeated  here  for  convenience 

X{k)    =  Y,    x(n)    e      "    >     0  ±  k  ±  N-l     .  (39) 

n=0 

The  product  formed  requires  four  real  multiplications  because 
both  the  exponential  and  the  input  data  are  complex  numbers, 
B  bits  in  length.  Each  multiplication  can  produce  a  result  of 
2B  bits  which  must  be  truncated  from  2B  to  B  bits,  hence  there 
are  four  quantization  errors  for  each  complex  valued 
multiplication.  The  four  quantization  errors  can  be  described 
as  four  independent  noise  sources.  The  truncation  noise  power 
is  defined  by  [Ref.  7]  as 
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a'T  =    [±2-2B   +  -2"2fl} 

3       4  (40) 

where  o\   =  truncation  noise  power 

B  -  number  of  bits  in  product  after  truncation     . 


D.   SCALING  NOISE 

The  FFT  as  an  algorithm  processes  a  data  vector  of  length 
N  by  successive  passes  at  the  vector.  During  each  pass,  the 
algorithm  performs  N/2  butterfly  operations.  Each  butterfly 
retrieves  two  complex  numbers  from  memory,  performs  the 
butterfly  computation  and  returns  two  complex  numbers  to  the 
same  memory  address.  The  complex  numbers  returned  to  memory  by 
the  butterfly  can  have  a  greater  magnitude  than  the  two 
complex  numbers  initially  fetched  from  memory  by  the 
butterfly.  In  order  for  the  results  of  the  butterfly  operation 
to  fit  in  the  fixed  word  length  of  the  memory  of  the  computer, 
the  results  must  be  truncated  (scaled)  .  Scaling  is  performed 
by  dividing  by  two  the  entire  data  array,  and  is  implemented 
in  software  by  a  right  shift  and  an  increment  of  the  arrays 
exponent  register. 

During  the  (m+1)  pass  of  the  data,  the  butterfly  algorithm 
selects  two  data  points  A(m)  and  B (m)  and  returns  to  memory 
A(m+1)  and  B(m+1).  The  truncation  error  results  from  the 
addition  of  two  numbers  of  like  sign  in  the  upper  part  of  the 
butterfly  algorithm  or  the  subtraction  of  two  numbers  of 
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opposite  sign  in  the  lower  part  of  the  butterfly  algorithm. 
The  input  to  the  butterfly  is  fixed  to  have  maximum  real  and 
imaginary  values  of  ±1 .  The  output  of  the  butterfly  is  fixed 
to  have  a  maximum  value  of  ±2,  double  the  input  value.  To 
prevent  overflow  of  the  butterfly  output,  either  a  prescale  by 
1/2  is  required  prior  to  entering  the  butterfly  or  an  extra 
bit  must  be  available  in  memory  to  accommodate  the  possible 
gain  of  two. 

One  method,  automatic  prescaling,  truncates  the  least 
significant  bit  of  the  butterfly  output  and  represents 
processing  noise.  In  the  case  for  which  no  scaling  is 
necessary  for  the  pass,  this  truncation  represents  significant 
processing  error  (noise) . 

A  second  method,  data  dependent  scaling,  is  performed  only 
if  a  butterfly  in  the  data  pass  overflows  without  scaling.  If 
any  word  in  memory  has  an  overflow  bit  set,  all  words  are 
shifted  right  as  they  are  fetched  from  memory  for  the  next 
pass  of  the  FFT.  The  total  number  of  right  shifts,  p,  executed 
during  the  FFT  process  is  used  at  the  output  of  the  FFT  as  a 
scale  factor  of  2P.  This  scale  factor  is  used  at  the  output  of 
the  FFT  so  that  the  total  processing  gain  is  maintained. 
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The  noise  power  created  by  the  scaling  of  the  data  is 
defined  by  [Ref.  7] 


ol 


—  {—2'2{B'1)    *   —2'2(B~1) } 


2      3  4 

where  a2s  =  scaling  noise  power 


(41) 


E.   OUTPUT  SIGNAL  TO  NOISE  RATIO 

Figure  21  shows  an  error  model  of  a  butterfly  in  an  FFT 


Figure  21.   FFT  error  model. 
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The  input  data  terms  to  a  butterfly  can  be  described  by 

Aim)    =  Aim)    +  NAim) 
Bim)    =  B{m)    +  NBim) 

(42) 

where  NA,   NB   =  additive  noise  terms 

Aim)  ,   Bim)    -  input  data  terms 

Aim)  ,   Bim)    -  noise  corrupted  data     . 

Without  the  effects  of  noise  the  butterfly  algorithm  produces 

two  outputs  described  by 

A(m+1)    =  Aim)    +  W(m)B{m) 

fl(/7?+l)  =  Aim)    -  Wim)Bim)  (43) 

where  Wim)    -  twiddle  factor     . 

The  error  due  to  the  A  component  will  be  calculated.  The  error 
due  to  the  B  component  is  statistically  equivalent  to  the 
error  in  the  A  component.  Including  the  effects  of  noise,  the 
output  of  the  butterfly  can  be  described  by 

Aim+1)    =  [Aim)    +  NAim)    +qsim)}    +  [Bim)    +  NBim)    +qs(m)]- 
[W{m)    +  qwim)  ]  +  qTim) 

where  qsim)    =  scale  noise 

qwim)    =  sine/ cosine  table  noise 

qTim)    =  truncation  noise   ( in  multiply) 

(44) 

The  total  power  at  the  output  of  the  butterfly  is  then 
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EiA(m+l)  A*(m+1)}   =  E[A(m)A*(m)  )    +  E{B(m)  B*  (m)  }    + 

E{NA{m)Nl{m)  }    +  E{NB{m)N*B{m)}    + 

(45) 
£{gs(/n)  gs*  (/n)  }    +  £{gs(m)  qfe  (m)  }    + 

E{S(in)S*(in)}E{grv(flj)gi;(/n)}    + 

E{qT{m)  qj  {m)  }      . 

Redefining   the   terms    of  the   total   power    leads    to 

E{A{m+l)  A' im  +  1)  }    =  S2  {m)    +  S2  im)    +  n\  (/n)    +  Nl  {m)    + 

2o\{m)    +  S2(m)ozw(m)    +  o|(/n) 

where  S2  im)  =  input  signal  power  '   ' 

o|(in)  =  scaling  noise  power 

o%(m)  -   truncation  noise  power 

ow(m)  =  sine/ cosine  table  noise  power     . 

Combining  like  terms  results  in 


£{A(/n+l)A*(/n+l)  }  =  2S2  (/n)  +  2iV^  (/n)  +205(17?)  + 

S2  (/n)  o^dn)  +  Or(/n)  . 


(47) 


This  equation  describes  the  total  power  as  a  function  of  the 
signal  power  and  the  individual  independent  uncorrelated  noise 
terms.  The  total  power  output  from  the  butterfly  can  also  be 
described  as  a  function  of  its  input  signal  power  and  the 
total  additive  noise  power  generated  within  the  butterfly  as 

E{A(m+l)A*  (m+1)  }    =  E{A(m+l)A*  (m+1)  }    +  (48) 

E{NA(m+l)Nx{m+l) }     . 

Redefining  these   terms 

E{A(m+l)A*  (m+1)  }    =S2(/n+l)    +2YJj(in+l).  <49) 
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Equation  47  and  4  9  are  equivalent.  The  signal  and  noise  power 
for  the  (m+1)  pass  can  be  formed  from  those  terms  in  the 
proceeding  pass 


S2  (m+1)    +  NJi(m+l)    =  2S2(m)    +  2n\  (m)    +  2o'Um)    + 

S2  (m)  o2w{m)    +  o2T{m)  . 


(50) 


Equating  the  signal  terms  in  Equation  50  yields 

S2(m+1)    =  2S2(m)   .  (51) 

Equating  the  noise  terms  in  Equation  50  yields 

NJi{m+l)    =  2Nl(m)    +  2a2s(m)    +  S2  (m)  o2w(m)    +  o2T(m)   .    (52) 
Equation  51  can  be  rewritten  to  show  how  the  input  signal 
power  increases  as  a  function  of  m  as  it  passes  through  the  mth 
stage  of  a  butterfly  of  an  FFT 

S2(m)    =  2mS2, 

(53) 

here  we  modeled  S2  (0)  =  S2 . 
The  signal  power  input  to  the  butterfly  is  defined  as 

S2   =  ±2-2{B-b)      .  (54) 

9 

Rewriting  Equation  52  to  include  this  exponential  increase  of 
S2  (m)  gives 
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NJlim+1)  =  2iVJ(/n)  +  2a2s(m)  +  2mS2o2w{m)  +  o2r(m)  .  (55) 
Equation  55  shows  the  relationship  between  the  noise  power  of 
a  butterfly  stage  as  a  function  of  the  noise  power  from  the 
proceeding  butterfly  stage.  The  noise  and  signal  power  can  now 
be  calculated  using  Equations  53  and  55  for  m  passes  through 
a  butterfly.  Once  calculated,  the  signal  to  noise  ratio  can  be 
formed. 

F.   RIGHT  JUSTIFIED  DATA 

The  length  of  a  word  in  the  computer  is  defined  to  be  16 
bits  (B)  and  accepts  a  12  bit  (b)  data  word.  Right  justified 
data  places  the  least  significant  bit  (LSB)  of  the  b  bit  data 
word  in  the  LSB  of  the  B  bit  processor  word.  The  output  noise 
to  signal  ratio  is  now  calculated  for  ten  passes  through  a 
butterfly  using  right  justified  data.  The  following 
assumptions  are  used  in  the  calculation; 

1.  The  input  signal  is  incoherent. 

2.  Scaling  is  performed  every  other  pass  after  the  most 
significant  bit  (MSB)  of  data  reaches  the  next  to  MSB 
of  the  processor. 

3.  The  first  two  passes  do  not  use  the  trigonometric 
table. 

After  the  tenth  pass  (i.e.,  FFT  size  of  1024),  it  is  shown 
(Appendix  D)  that  the  noise  to  signal  power  is 
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B  and  b  are  previously  assumed  to  be  16  and  12  bits  in  length 

respectively.  Substituting  these  values  into  Equation  56 


(57) 


The  noise  to  signal  ratio  reduces  to 

AT2 

-^(10)  =  -67.5  db     -   95.1  dB     -71.1  dB     -79 . 1  dB  (58) 

=  input  trig  truncation      scaling 

noise  noise  noise  noise     . 

Clearly,  the  dominant  degradation  noise  is  the  truncation 
table  noise  (-71.1dB),  followed  by  the  scaling  noise  to  signal 
ratio  (-79.1dB)  .  For  right  justified  data,  the  signal  to  noise 
ratio  is  calculated  to  be  65.71dB. 

G.   LEFT  JUSTIFIED  DATA 

Left  justified  data  places  the  most  significant  bit  (MSB) 
of  the  b  bit  data  word  in  the  MSB  of  the  B  bit  FFT  processor 
word.  The  output  noise  to  signal  ratio  is  now  calculated  for 
10  passes  through  a  butterfly  for  left  justified  data.  The 
following  assumptions  are  made 
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1.  The  signal  is  incoherent. 

2.  Data  is  left  justified  on  input  and  scaling  is 
performed  every  other  pass  starting  with  the  first 
pass . 

3.  The  first  two  data  passes  do  not  use  the 
trigonometric  table. 

It  can  be  shown  that  after  the  tenth  pass  through  the 
butterfly  that  the  noise  to  signal  ratio  is 

^(10)  =  ^  +  4o2„  +  1^1  ♦  2-^.  (59) 

S2  S2  4  s2  S2 

B  and  b  were  previously  assumed  to  be  16  and  12  bits 
respectively.  Substituting  these  values  into  Equation  59  we 
obtain 


—  (10) 
S2 


The  noise  to  signal  ratio  reduces  to 


N2 

—  (10)  =  -67.5  dB     -95.1  dB     -89 . 1  dB     -77  . 1  dB  (61) 

S2 

=  input  trig  truncation       scaling 

noise  noise  noise  noise  . 


Clearly,  the  worst  noise  contribution  to  the  SNR  is  the 
scaling  noise  (-77.1dB),  followed  by  the  truncation  noise  (- 
89.1dB) .  For  left  justified  data,  the  signal  to  noise  ratio  is 
calculated  to  be  67.0dB. 
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Further  analysis  shows  that  for  large  FFT's  (i.e.,  in  the 
range  of  215  and  above) ,  the  trigonometric  table  noise  is  the 
dominant  noise  source.  A  simulation  was  performed  using  an 
AMDAHL  machine  to  examine  the  error  generated  in  performing  a 
FFT  and  an  iFFT.  Figure  22  shows  the  experimental  setup. 


Data 

AMDAHL  Machine 

Data 

FFT 

».       i  CCT 

*      ihr  1 

x\ 

save  error  in  table 

Figure  22.  FFT  processing  error  flowchart  [Ref.  8]. 

The   error  that  was   generated  was  tabulated.   Figure  23 
illustrates  the  error  as  a  function  of  FFT  length. 


B.   SUMMARY 

For  either  left  or  right  justified  data,  the  dominant 
noise  at  the  output  of  the  1024  point  FFT  is  the  scaling  or 
truncation  noise.  The  SNR  at  the  output  of  the  FFT,  for  left 
justified  data  is  67.0dB  and  exceeds  the  SNR  for  right 
justified  data  which  is  65.71dB.  Clearly,  left  justified  data 
maximizes  the  computational  SNR  at  the  output  of  the  FFT  for 
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Figure  23 
[Ref.  8]. 


FFT  error  as  a  function  of  transform  length 
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the  values  of  b,  B  and  N  used  in  the  above  examples. 

For  a  right  justified  data  set  transform  size  of  N=213 
(8192  points) ,  the  noise  to  signal  ratio  after  the  13th  pass 
can  be  written 

Nl  213N*   +    (215  +213  +  212)S2a2w   +  2047<4  +  1700^   (62) 

I2  21352 

This  equation  reduces  to 

—  13   =  — i  +  5.5ol   +  —  — -  +  — = .       <63> 

S2  S2  4  s2  48.2  52 

Comparing  Equation  5  6  to  Equation  63,  we  see  that  by 
increasing  the  transform  size  from  210  to  213  increases  the 
noise  to  signal  ratio.  The  trigonometric  table  noise  is 
increasing  faster  than  either  the  scaling  or  truncation  noise. 
The  signal  to  noise  ratio  is  reduced,  from  65.71dB  to  65.68dB. 
As  the  transform  size  is  increased  to  accommodate  larger  data 
sets,  the  trigonometric  table  noise  becomes  the  dominant  error 
term. 
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VII.  SUMMARY 

A.   CONCLUSIONS 

A  frequency  domain  based  algorithm  was  developed  and 
tested  to  estimate  the  differential  arrival  time  of  a  pulsed 
radar  signal  collected  by  two  passive  sensors.  For 
convenience,  the  actual  implementation  was  performed  through 
time  domain  processing  even  though  frequency  domain  processing 
is  advocated.  The  performance  of  the  algorithm  is 
characterized  by  a  ROC  curve  as  a  function  of  SNR.  For  3dB  SNR 
and  a  Pf(i  of  0.01,  the  probability  of  making  a  correct  TDOA 
estimate  exceeds  that  of  an  incorrect  estimate  for  a  single 
look.  At  18dB  SNR,  the  probability  of  making  a  correct  TDOA 
estimate  is  100  percent  for  all  P£a's.  If  multiple  looks  (i.e., 
multiple  bursts)  are  allowed,  the  probability  of  making  a 
correct  TDOA  decision  at  each  SNR  will  increase. 

An  I/Q  demodulator  is  assumed  in  the  radar  receiver.  The 
pdf  of  the  signal  driving  the  correlation  detector  is 
determined  from  where  it  is  obtained  in  the  I/Q  receiver.  The 
pdf,  depending  on  selection,  will  be  zero  mean  Gaussian  versus 
non-zero  mean  Gaussian,  or  chi-squared  versus  non-central  chi- 
squared,  or  Rayleigh  versus  Rician  (non-central  Rayleigh) . 
This  thesis  assumes  an  envelope  detector  at  the  output  of  the 
I/Q  demodulator.  The  envelope  detector  output  has  a  Rayleigh 
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or  Rician  pdf  and  drives  the  correlation  detector.  The 
calculation  of  a  CFAR  threshold  assumes  a  Gaussian  pdf  at  the 
output  of  the  correlation  detector.  For  small  time  lags, 
summation  of  the  output  of  the  correlator  may  produce  a 
probability  distribution  that  deviates  from  Gaussian.  This 
deviation  would  produce  a  biased  threshold  calculation.  To 
minimize  the  bias,  a  sufficient  number  of  terms  must  be  summed 
at  the  output  of  the  correlator  to  allow  the  Gaussian 
approximation . 

If  the  received  pulses  are  collected  in  the  frequency 
domain,  a  spectral  domain  based  correlation  algorithm  can  be 
implemented.  An  algorithm  is  given  in  this  thesis. 

For  any  digital  signal  processing  algorithm  that  uses  the 
FFT,  processing  errors  must  be  considered.  For  large  transform 
sizes  the  trigonometric  noise  power  dominates.  The  length  of 
the  trigonometric  coefficient  word  affects  the  degradation  of 
the  SNR  at  the  output  of  the  FFT.  The  larger  this  word  the 
smaller  the  noise  power.  For  a  given  transform  size,  left 
justified  data  will  have  a  higher  output  SNR  than  right 
justified  data. 

For  a  correct  TDOA  estimate,  the  correlation  algorithm 
requires  the  reception  of  the  leading  pulses  in  the  radar 
pulse  burst.  If  the  pulses  are  received  in  an  adequate  SNR, 
the  correlation  algorithm  will  produce  a  well  defined  peak 
allowing  the  TDOA  estimate.  Should  the  environment  degrade  the 
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received  signal  (i.e.,  destructive  interference  due  to 
multipathing,  weather  or  terrain) ,  the  correlation  algorithm 
should  be  discarded  in  favor  of  the  traditional  angle  of 
arrival  (AOA)  algorithm  (i.e.,  maximize  received  energy) . 

If  the  radar  uses  a  staggered  PRF,  the  reception  of  the 
first  pulse  or  several  pulses  has  little  impact  on  the 
position  of  the  correlation  peak.  Under  this  condition,  the 
TDOA  algorithm  is  superior  to  the  AOA  algorithm. 

B.   RECOMMENDATIONS 

A  ROC  curve  should  be  constructed  for  the  correlation 
algorithm  used  in  this  thesis  for  multiple  TDOA  looks.  An 
intensity  display  could  be  designed  that  would  plot  as  a 
function  of  time  (i.e.,  snapshots),  those  lag  points  chosen  as 
TDOA  estimates.  Patterns  displayed  on  the  intensity  plots 
would  allow  the  user  to  visually  determine  the  correct  TDOA 
(i.e.,  incoherent  averaging). 

The  performance  of  the  spectral  domain  based  correlation 
algorithm  developed  in  this  thesis  must  be  further  quantified. 
Sets  of  ROC  curves  should  be  obtained  for  different  SNR' s 
(i.e.,  SNR  variations  between  channels). 

The  simulated  FFT  error  graph  should  be  validated  by 
measuring  the  error  performance  of  a  dedicated  properly 
dimensioned  FFT. 
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APPENDIX  A.  CORRELATION  MEAN  AND  VARIANCE 

Two  zero  mean,  independent  time  series  are  the  inputs  to 
a  correlator.  This  appendix  will  derive  the  expected  value  and 
variance  of  the  crosscorrelation  function  output.  In  this 
thesis,  the  noise  is  zero  mean  (i.e.,  shifted)  Rayleigh  noise. 

The  crosscorrelation  function  rIy(l)  is  defined  in  Equation 
8  and  is  repeated  here  for  convenience 


N-l-\l 

E 

i=0 


r*y(i)  =   £  x(i)y(i   +  1)      .  (64) 


A.   EXPECTED  VALUE 

x±  and  y±  are  assumed  to  be  independent,   zero  mean 
sequences 

E{x±)    =  E{y±)    =  0.  (65) 

The  expected  value  of  rxy(l)  is 

W-1-|I|  JV-1-|1| 

Eir^d)}    =  E{    £     x,yi+i}    =     £     ^Cxjyj+J} 

i=0  i=0 

w-iHil  (66) 

1=0 

=  0      . 
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B .   VARIANCE 

The  random  variable  z  has  a  variance  defined  to  be 

o\  =  E{{z   -  Eiz})2} 
because  E{z}    =  0  for  a  zero  mean  sequence,  (<>7) 

o2z  =  E{z2}  . 

The  variance  of  the  crosscorrelation  function  rxy(l)  is  then 

w-i-UI  w-i-|i|      w-i-|i| 

o2rxyU)  =  E{r2xy)  *  eh  £  *i^i>2>  -  *<  E  x±y±*i  E  *&±+$ 

i=o  I=o        i=o       (6g) 

w-i-|i|  «r-i-|i| 

=  E{  E    E  x^i*i  xjyj+i]    • 

1*0    j=0 

x^  and  yi+1  yj+1  are  two  groups  that  are  independent  of  each 
other.  The  expectation  operation  can  be  applied  to  each  group 
individually . 

V-l-|l|  W-1-|I| 

<yu>  =  E    E  BiXjXj)   Eiy^y^}      .  (69) 

1=0      J=0 

The  two  indices  i  and  j  define  a  matrix  of  values.  Two 
contributions  have  to  be  considered. 

Contribution  1.  When  i=j  the  summation  occurs  along  the 
diagonal  of  the  matrix.  Only  one  summation  is  used  and 
essentially  the  terms  are  being  squared.  These  squared  terms 
are  not  independent.  To  simplify  the  computation  define  m 
equal  to  i  and  j 
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q-"£       E{x'}E{y2+1}      .  (70) 


m=o 


(71) 


Because    x„   and   yn   are    zero   mean    sequences, 

E{x2}    =  ol,    Eiy^}    =  a2y 

and  o\  -  o2y  =  o2     . 
Substituting   a2    into   Equation   70 

N-l-\l\  W-1-|I| 

q  =     £    oioj  =  o<   £    i 

m=0  m=0 

=   o4[  (iV-l-|i|)  +1] 
=   04(W-|I|)       . 

Contribution  2.  When  i*j  all  the  terms  except  those  on  the 
diagonal  are  summed.  These  terms  are  independent.  For  i*j 

tf-i-UI  w-i-|i| 
c2  =  E    E  EixjEixjEiy^jEiy^j) 

i=o    £o  (72) 

=  0,  using  Equation  6  5 


Therefore,  Equation  69  becomes  a2.   =  C,  +  C, 

*  "<"  1         2  (73) 

=  a4(N  -  1 1 1  )  . 
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C.   CROSS CORRELATION  STATISTICS 

Each  data  point  output  by  the  correlation  function 
represents  a  summation  of  terms.  Using  the  central  limit 
theorem,  the  output  of  the  correlation  function  is  assumed  to 
have  a  Gaussian  distribution. 


In  conclusion  r^d)    ~N(0,   a4  ( AT- 1 J  I ) 

where  N  =  Normal    {Gaussian)   distribution     . 
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APPENDIX    B.     MATLAB    CODE 

A.  USERS  GUIDE 

The  frequency  domain  based  correlation  software  given  in 
this  appendix  is  written  in  Pro-Matlab  (version  3.5h).  All 
simulations  were  run  on  a  Naval  Postgraduate  School  Sun  work 
station  using  a  UNIX  operating  system. 

B.  CORRELATION  MATLAB  CODE 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%    Program  Name  :  FreqCorr7.m  12  May  91 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

clg 

clc 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  PULSE  CORRELATION 

%  Version  7.0 

% 

%  1 .  1  Channel  vs  Reference  Channel 

%  or 

%     1  Channel  vs  Second  Channel . 

% 

%  2.  Rician  Signal  and  Rayleigh  Noise  pdf ' s  . 

%     (Low  SNR  signals  have  power  approximately  2.) 

% 

%  3.  0  dc  Correlation.  Subtract  dc  from  both  signals. 

% 

%  4.  Normalize  correlation  output  in  Fourier  domain. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  This  MatLab  program  will  either  : 

%  1.  Correlate  a  received  pulse  sequence  against  a  reference 

%   pulse  sequence.  The  reference  sequence  parameters  are 

%    specified  by  the  user  and  are  used  to  construct  the 

%    sequence. 

%  or 
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%  2.  Correlate  two  received  pulse  sequences  against  each 
%     other. 


%   Correlation  of  the  pulse  sequences  is  accomplished  by 
%   multipling  in  the  frequency  domain  one  spectrum  against 
%   the  conjugate  of  the  other  spectrum.  An  inverse  FFT  is 
%   performed  on  the  result  to  yield  the  time  difference 
%   (TDOA)  between  the  two  sequences. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clc 

echo  on 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  SELECT  CORRELATION  OPTION 

% 

%     Select  either  1  or  2  : 

% 

%   1.  Time  delay  in  one  receiver  channel  measured  against  a 

%      reference  signal. 

% 

%   2.  Differential  time  delay  between  two  channels. 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

echo  off 

option  =  input  (' Select  correlation  option  1  or  2  ;       '); 

clc 


echo  on 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Option  1  :  Time  Delay  In  One  Channel 

% 

%         Reference  Pulse  Train  Parameter  Selection. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

echo  off 

if  option==l 

wi=input (  'Pulse  Width  (seconds)  :     '); 

pr=input (  'Pulse  Period  (T  seconds)  :     '); 

nu=input  (  'Number  of  pulses  in  pulse  train  :     '); 

NumberPointsRef=pr*nu; 
end 
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clc 

echo  on 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% 

%         Delayed  Pulse  Train  Parameter  Selection 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

echo  off 

dly=input (' Channel  1  delay  tau  (sec.)  of  pulse  train  :    '); 

wi2=input  (' Pulse  Width  (seconds)  :    '); 

pr2=input (' Pulse  Period  (T  seconds)   :    '); 

nu2=input (' Number  of  pulses  in  Delayed  Pulse  Train  :    ')/ 

NumberPointsDel= (pr2*nu2) +dly;   %calc.  nmbr  pts  sig  +  delay 

clc 

if  option==2 

dlyCh2=input (' Channel  2  delay  tau  (sec.)  pulse  train  :    ')/ 

wiCh2=input (' Pulse  Width  (seconds)  :    '); 

prCh2=input ('Pulse  Period  (T  seconds)   :    '); 

nuCh2=input (' Number  of  pulses  in  Delayed  Pulse  Train  :    '); 

NumberPointsDel2= (prCh2*nuCh2) +dlyCh2;  %nmbr  pts  sig  +  delay 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   Zero  pad  both  sequences  to  N3=N2+N1-1.  Must  then  make  N3 
%   meet  the  simple  formula  N3  =  2Am  to  allow  speedy 
%   computation  of  the  FFT. 

%   WARNING  :  IBM  80286  MATLAB  WILL  NOT  ALLOW  VECTORS  GREATER 

%   THAN  4000.  So  an  input  pulse  train  of  128  will  exceed  the 

%   system.  Solution  is  to  use  the  SUN  workstation  or  386/486. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

if  option==l 

ZeroPadPoints=NumberPointsRef+NumberPointsDel-l; %N3=N2+N1-1 

for  m=2:l:19;  %2~  (19)  =  524288  point  FFT  max 

if  2*m  >=  ZeroPadPointsf 

ZeroPadPoints=2Am; 

break;      %when  find  a  2/sm  power,  exit  loop 
end 
end 
end 
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if  option==2 

ZeroPadPoints=NumberPointsDel2+NumberPointsDel-l; %N3=N2+N1-1 
for  m=2:l:19;  %2"  (19)  =  524288  point  FFT  max 

if  2~m  >=  ZeroPadPoints, 

ZeroPadPoint s=2Am; 

break;        %when  find  a  2^  power,  exit  loop 
end 
end 
end 


%%%%%%%%%%%%%%%  Create  Reference  Pulse  Train   %%%%%%%%%%%%%% 
if  option==l 

ReferencePulse=zeros (1 : ZeroPadPoints) ;  %zero  pad  past  pr*nu 
for  PulseCounter=l :pr :NumberPointsRef   %steps  of  Ref.  period 
for  CutUpPulse=0 :pr; 
if  CutUpPulse<=wi 

ReferencePulse (PulseCounter+CutUpPulse) =  1.0; 
end 
end 
end 

Ref dc=mean (ReferencePulse (1 :NumberPointsRef) ) ; 
ReferencePulse (1 :NumberPointsRef ) =ReferencePulse (1 :NumberPoi 
ntsRef ) -Refdc; 
end 


%%%%%%%  Create  Channel  1  Delayed  Pulse  Train   %%%%%%%%%%%%%% 

DelayedPulseTrain=zeros (1 : ZeroPadPoints) ; %0  pad  past  pr2*nu2 

for  PulseCounter= (dly+1) :pr2: (pr2*nu2+ (dly+1) ); %pr2*nu2=Nmbr 

for  CutUpPulse=0 :pr2;         %...pts  in  Del.  Pulse  Train 

if  CutUpPulse<=wi2 

DelayedPulseTrain (PulseCounter+CutUpPulse) =1.0; 
end 
end 
end 


if  option==2 

%%%%%%%  Create  Channel  2  Delayed  Pulse  Train   %%%%%%%%%%%%%% 
DelayedPulseTrain2=zeros (1 : ZeroPadPoints) ; %0  padpast  pr2*nu2 
for  PulseCounter=(dlyCh2+l) :prCh2: (prCh2*nuCh2+ (dlyCh2+l) ) ; 
for  CutUpPulse=0:prCh2; 
if  CutUpPulse<=wiCh2 

DelayedPulseTrain2 (PulseCounter+CutUpPulse) =1.0; 

74 


end 
end 
end 
end 


clc 

echo  on 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Signal  to  Noise  Ratio  Calculation 

%  Define  Rayleigh  noise  power  to  equal  2. 

%       Rician  pdf  created  for  sinusoid  +  Gaussian  noise. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

echo  off 

SNR=input('  Desired  SNR  :   ')/ 

NoisePwr=2;  %rayleigh  noise  pwr  =2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%       Sum  Channel  1  &  2  signals  plus  their  initial  delay. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SigEnergy=0;  %initialize  Channel  1  energy  to  0 

for  index=l :NumberPointsDel; %sum  over  Ch.  1  delayed  signal 

SigEnergy=SigEnergy+ ( (DelayedPulseTrain (index) ) . ~2) ; 
end 

if  option==2 

SigEnergy2=0;  %initialize  Channel  2  energy  to  0 

for  index=l :NumberPointsDel2; %sum  over  Ch.  2  delayed  signal 
SigEnergy2=SigEnergy2+ ( (DelayedPulseTrain2 (index) )  . A2)  ; 
end 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Calc.  Ch.  1  signal  amplitude  to  meet  required  input  SNR. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SNR10=SNR/10; 

SNRLinear=10ASNR10;  %SNR  in  linear  units 

SigAmplitude=sqrt (NoisePwr*SNRLinear*NumberPointsDel/SigEner 

gy)  / 


if  option==2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Calc.  Ch.2  signal  amplitude  to  meet  required  input  SNR. 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SNR10=SNR/10; 

SNRLinear=10~SNR10;  %SNR  in  linear  units 

SigAmplitude2=sqrt (NoisePwr*SNRLinear*NumberPointsDel2/SigEn 

ergy2) ; 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Create  Rayleigh  Noise  matrix  for  Channel  1. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
rand (' normal' ) ;  %Gaussian  raean=0,  variance=l 

for  index=l :NumberPointsDel, 

GaussNoisel (index) =rand;      %N (0, 1) 

GaussNoise2 (index) =rand;      %N  (0, 1) 
end 
RayleighNoise=sqrt ( (GaussNoisel) . ^2+ (GaussNoise2) . ^2) / 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%   Create  Rayleigh  Noise  matrix  for  Channel  2 . 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

if  option==2 

rand (' normal' ) ;  %Gaussian  mean=0,  variance=l 

for  index=l : NumberPointsDel2 , 

GaussNoise3 (index) =rand;      %N (0, 1) 

GaussNoise4 (index) =rand;      %N (0, 1) 
end 

RayleighNoise2=sqrt ( (GaussNoise3)  . A2+ (GaussNoise4)  . A2)  / 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   Make  Rician  amplitude  of  Channel  1  delayed  pulse  train 
%   and  remove  the  dc  component  of  entire  signal. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SE=0; 

for  index=l :NumberPointsDel,   %1  -->  delay  +  delayed  signal 
if  DelayedPulseTrain (index) ==1        %if  one,  make  Rician 


DelayedPulseTrain (index) =DelayedPulseTrain (index) . *rician (Si 
gAmplitude) / 

SE= (DelayedPulseTrain (index) ) . A2+SE; 
end    %end  Rician  modification  loop 
end 
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SP=SE/NumberPointsDel; 

%%%%%%%   Add  noise  to  Channel  1  Delayed  Puxse  Train  %%%%%%% 
for  index=l :NumberPointsDel,   %1  -->  delay  +  delayed  signal 
if  DelayedPulseTrain (index) ==0  %no  Rayleigh  noise  to  signal 

DelayedPulseTrain (index) =DelayedPulseTrain (index) +RayleighNo 
ise (index) ; 

end 
end 

%%%%%  Remove  dc  component  of  Channel  1  Pulse  Train  %%%% 
Chldc=mean (DelayedPulseTrain (1 :NumberPointsDel) ) / 
DelayedPulseTrain (1 :NumberPointsDel) =DelayedPulseTrain (1 :Num 
berPointsDel) -Chide; 


if  option==2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Make  Rician  amplitude  of  Channel  2  delayed  pulse  train 
%  and  remove  dc  component  of  entire  signal. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SE=0; 

for  index=l :NumberPointsDel2,  %1  -->  delay  +  delayed  signal 
if  DelayedPulseTrain2 (index) ==1        %if  one,  make  Rician 

DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) . *rician ( 
SigAmplitude2) ; 

SE=  (DelayedPulseTrain2  (index)  )  .  /S2+SE; 
end   %end  Rician  modification  loop 
end 
SP=SE/NumberPointsDel2  ; 


%%%%%%   Add  noise  to  Channel  2  Delayed  Pulse  Train  %%%%%%%% 

for  index=l :NumberPointsDel2,  %1  -->  delay  +  delayed  signal 

if  DelayedPulseTrain2 (index) ==0  %no  Rayleigh  noise  to  sig 

DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) +Rayleigh 
Noise2 (index) / 

end 
end 

%%%%%  Remove  dc  component  of  Channel  2  Pulse  Train  %%%% 
Ch2dc=mean (DelayedPulseTrain2 (1 :NumberPointsDel2) ) / 
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DelayedPulseTrain2 (1 :NumberPointsDel2) =DelayedPulseTrain2 (1 : 

NumberPointsDel2) -Ch2dc; 

end 

subplot (221)  %two  rows,  two  columns 

time= (0 : 1 : ZeroPadPoints-1) /      %Start  time  axis  at  0 

if  option==l 

%%%%%%%%%%%%%  Plot  Reference  Pulse  Train   %%%%%%%%%%%%%%%%%% 

plot (time, ReferencePulse) ; 

title (' Reference  Pulse  Sequence' ) ; 

xlabel ( ' time' ) / 

grid 

end 


%%%%  Plot  Channel  1  pulse  train  +  Rayleigh/Rician  noise%%%% 

topDel=l . l*max (DelayedPulseTrain) ; 

bottomDel=l . l*min (DelayedPulseTrain) ; 

axis([0  NumberPointsDel  bottomDel  topDel] ) 

plot (time, DelayedPulseTrain) ; 

title (' Channel  1  Pulse  Sequence  +  Rician/Rayleigh  Noise')/ 

xlabel (' time' ) ; 

grid 


if  option==2 

%%%%  Plot  Channel  2  pulse  train  +  Rayleigh/Rician  noise  %%% 

topDel2=l . l*max (DelayedPulseTrain2) ;  %10  percent  headroom 

bottomDel2=l . l*min (DelayedPulseTrain2) ; 

axis([0  NumberPointsDel2  bottomDel2  topDel2]) 

plot (time, DelayedPulseTrain2)  ; 

title (' Channel  2  Pulse  Sequence  +  Rician/Rayleigh  Noise'); 

xlabel ( ' time' ) / 

grid 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Correlation 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

if  option==l 

fR=f ft (ReferencePulse) ; 

fD=f ft (DelayedPulseTrain) / 

end 
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if  option==2 

fR=f ft (DelayedPulseTrain) / 

fD=fft (DelayedPulseTrain2)  ; 

end 

z=fD. *conj (fR) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Calculate  normalizing  coefficient 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Norml=0; 

Norm2=0 

for  k=l:l:length(fD) 

Norml= (abs (fD (k) )  .  ~2) +Norml; 

Norm2=(abs (fR(k) ) . A2) +Norm2; 
end 

Norml=sqrt (Norml) ; 
Norm2=sqrt (Norm2) ; 

NormCoef f=Norml*Norm2; 
z=z . /NormCoef f; 

z=z* length (fD) ; 

ifftz=ifft (z)  ; 
magz=abs (ifftz)  ; 


%%%%%  Flip  vector  for  appearance  %%%%%%%%%%%%%%% 

Topmagz= (length (magz) /2+1) ; 

Bottommagz=length (magz) /2; 

magzLen=length (magz) ;   %measure  length  of  magz 

Transmagz= [magz (Topmagz :magzLen) ,magz (1 :Bottommagz) ] ; 

%%%%%%  Make  -time  to  +time  axis  %%%%%%%%%%%%%%%% 
minustime=- (length (Transmagz) /2) / 
posittime= (length (Transmagz) /2) -1; 
timel= (minustime: 1 :posittime) ; 

topmagz=l. l*max (magz) ;  %add  10  percent  head  room 
axis ( [minustime  posittime  0  topmagz])  %scale  axis 

subplot  (212) 

plot (timel, Transmagz) ; 

title('FFT  Correlation  :  Channel  1  vs.  Channel  2'); 
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xlabel ( ' time  lag')/ 

grid 

gtext (' +7dB  SNR' ) 

lj4print  %Spanagel  Room  427  Laserjet 

axis([l  2  3  4]);   %reset  axis  scaling 
axis; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Program  Name  :  CorrPfaLoop.m  19  July  91 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Purpose  :  To  compute  the  Pfa  for  Signal  /  Noise  only  for 
%  then  correlation  function.  This  program  creates  zero  mean 
%  RAYLEIGH  noise. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

clg 

subplot (221) 

rand ( ' normal'  ) 

PFA=1.0;  %will  be  divided  by  ten  to  start 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for  PFAC0UNTER=1 : 1 : 4 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

PFA=PFA/10;  %Pfa  =  0.1  0.01,  0.001,  0.001 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for  SNR=0:1:18  %Outer  loops  through  all  SNR' s . 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Iterations=30;  %30  data  points  per  SNR  point 

%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for  loop=l : 1 : Iterations  %inner  loop,  creates  groups  of  +-  a 

%%%%%%%%%%%%%%%%%%%%%%%%%%% 

loop; 

AboveThreshold (loop) =0; 

BelowThreshold(loop) =0; 
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Figure  24.   FreqCorr7.m  MATLAB  flowchart 
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Threshold (loop) =0 ; 
MaxRxy ( loop) =0 / 
RxyThreelag  (loop) =0; 

%%%%%%%%%%%%  Set  Parameters  %%%%%%%%%%%%%%%%%%%%%%%% 
N=100;   %100  noise  points,  therefore  200  correlation  points 
a=10;    %+-  10  %  on  either  side  of  0  time  lag 
Pfa=PFA;     %Reestablish  because  Pfa  destroyed  each  loop 
PfaDef ined=Pfa; 

%%%%%%%%%%%%%  Make  N  points  Noise  %%%%%%%%%%%%%%%%%% 
snrnoiseloop  %call  signal  plus  noise  generation  at  a  SNR 


%%%%%  Correlation  of  zero  mean  noise  only  for  Pfa  calcul.  %% 
MeanRayl=mean (RayleighNoise) /    %noise  created  in  snrnoise.m 
MeanRay2=mean (RayleighNoise2) ; 
RayleighNoise=RayleighNoise-MeanRayl ; 
RayleighNoise2=RayleighNoise2-MeanRay2; 

Rxy=xcorr (RayleighNoise, RayleighNoise2) ;  %xcorr  signal+noise 
later 

%%%%%%%%%%%%%%%%%%  Compute  t~  %%%%%%%%%%%%%%%%%%%%%%% 

t=0; 

for  i=(N-a) :1: (N+a) 

t=t+Rxy (i)  ; 
end 
t=t/(2*a+l) ; 

%%%%%%%%%%%%%%%  Compute  variance  t  %%%%%%%%%%%%%%% 

Vart=0; 

for  i=(N-a) :1: (N+a) 

Vart=Vart+ (Rxy (i) -t)  ."2; 
end 
Vart=Vart/ (2*a+l) ; 

%%%%%%%%%%%%%%%%%  Define  Pfa  %%%%%%%%%%%%%%%%%%%% 

Pfa=Pfa*2;  %multiply  by  2 

Pfa=l-Pfa;  %subtract  one 

Threshold (loop) =inverf (Pfa) ; 

Threshold (loop) =Threshold (loop) *sqrt (2) *sqrt (Vart)  /  %mult  by 

the  sqrt(2) 

%%%%%  Xcorr  of  signal  +  zero  mean  noise  %%%%%%%%%% 
Rxy=xcorr (DelayedPulseTrain, DelayedPulseTrain2) ; 
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%%%%%%%  Compare  +a  to  -a  to  Threshold  %%%%%%%%%%%% 
for  i= (N-a) : 1: (N+a) 

if  Rxy (i) >=Threshold (loop) 

AboveThreshold (loop) =AboveThreshold (loop) +  1; 
else 

BelowThreshold (loop) =BelowThreshold (loop) +  1; 
end 

if  Rxy (i) >MaxRxy (loop) 

MaxRxy (loop) =Rxy (i) / 

SaveThreshold (loop) =Threshold (loop) ; 

%SaveVart (loop) =sqrt (Vart) ; 
end 
end 


RxyThreelag (loop) =Rxy (97) ; 


end 


%end  #  Iterations  loop 


%%%%%%%  Count  correct  threshold  crossings  %%%%%%%% 

Correct TDOA=0 ; 

for  count=l : 1 : Iterations 

if  RxyThreelag (count) ==MaxRxy (count) 
Correct TDOA=Correct TDOA+1 ; 

end 
end 
PercentCorrectTDOA= (Correct TDOA*100) /Iterations; 


diary  flag4 

PercentCorrectTDOA 

BelowThreshold 

AboveThreshold 

PfaDefined 

SNR 

SaveThreshold 

lag  below 

RxyThreelag 

three  pulses. 

MaxRxy 

diary  off 


end 
end 


%compare  threshold  to  zeroth 
%Channel  2  lags  Channel  1  by 

%end  #  test  points  loop 
%end  Pfa  loop 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%    Program  Name  :  snrnoiseloop .m  19  JULY  91 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

dly=0;    ^Channel  1  delay 

wi2=5;    %Channel  pulse  width 

pr2=10;   %Channel  1  pulse  period 

nu2=10;   %Channel  1  number  of  pulses 

NumberPointsDel= (pr2*nu2) +dly;   %calc.  nmbr  pts  sig  +  delay 

dlyCh2=3;   %Channel  2  delay 

wiCh2=5;    %Channel  2  pulse  width 

prCh2=10;   %Channel  2  pulse  period 

nuCh2=10;   %Channel  2  number  of  pulses 

NumberPointsDel2= (prCh2*nuCh2) +dlyCh2;  %nmbr  pts  sig  +  delay 

ZeroPadPoints=NumberPointsDel2+NumberPointsDel-l;%N3=N2+Nl-l 
for  m=2:l:19;  %2~(19)  =  524288  point  FFT  max 

if  2Am  >=  ZeroPadPoints, 

ZeroPadPoints=2'Nm/ 

break;        %when  find  a  2"m  power,  exit  loop 
end 
end 

%%%%  Create  Channel  1  Delayed  Pulse  Train   %%%%%%%%%%% 
DelayedPulseTrain=zeros (1 : ZeroPadPoints) ;%0  pad  past  pr2*nu2 
for  PulseCounter=(dly+l) :pr2: (pr2*nu2+ (dly+1) ) ; %pr2*nu2=Nmbr 
for  CutUpPulse=0 :pr2;         %...pts  in  Del.  Pulse  Train 
if  CutUpPulse<=wi2 

DelayedPulseTrain (PulseCounter+CutUpPulse) =  1.0; 
end 
end 
end 

%%%%%  Create  Channel  2  Delayed  Pulse  Train   %%%%%%%%%% 
DelayedPulseTrain2=zeros (1 : ZeroPadPoints) ; %0  padpast  pr2*nu2 
for  PulseCounter=(dlyCh2+l) :prCh2: (prCh2*nuCh2+ (dlyCh2+l) ) ; 
for  CutUpPulse=0 :prCh2; 
if  CutUpPulse<=wiCh2 

DelayedPulseTrain2 (PulseCounter+CutUpPulse) =  1.0; 
end 
end 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Signal  to  Noise  Ratio  Calculation 

%  Define  Rayleigh  noise  power  to  equal  2. 
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%        Rician  pdf  created  for  sinusoid  +  Gaussian  noise. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

echo  off 

NoisePwr=2;  %rayleigh  noise  pwr  =2 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%       Sum  Channel  1  &  2  signals  plus  their  initial  delay. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SigEnergy=0;  %initialize  Channel  1  energy  to  0 

for  index=l :NumberPointsDel; %sum  over  Ch.  1  delayed  signal 

SigEnergy=SigEnergy+ ( (DelayedPulseTrain (index) ) . "2) ; 
end 

SigEnergy2=0;  %initialize  Channel  2  energy  to  0 

for  index=l :NumberPointsDel2; %sum  over  Ch.  2  delayed  signal 
SigEnergy2=SigEnergy2+ ( (DelayedPulseTrain2 (index) ) . A2) / 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Calc.  Ch.  1  signal  amplitude  to  meet  required  input  SNR. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SNR10=SNR/10; 

SNRLinear=lCrSNR10;  %SNR  in  linear  units 

SigAmplitude=sqrt (NoisePwr*SNRLinear*NumberPointsDel/SigEner 

gy)  ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%   Calc.    Ch.2    signal   amplitude   to  meet    required   input    SNR. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

SNR10=SNR/10; 

SNRLinear=10,NSNR10;  %SNR  in   linear  units 

SigAmplitude2=sqrt (NoisePwr*SNRLinear*NumberPointsDel2/SigEn 

ergy2) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%    Create  Rayleigh  Noise  matrix  for  Channel  1. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for  index=l :NumberPointsDel, 

GaussNoisel (index) =rand;      %N (0, 1) 

GaussNoise2 (index) =rand;      %N (0, 1) 
end 
RayleighNoise=sqrt ( (GaussNoisel) . A2+ (GaussNoise2) . A2) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   Create  Rayleigh  Noise  matrix  for  Channel  2. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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for  index=l :NumberPointsDel2, 

GaussNoise3 (index) =rand;      %N (0, 1) 

GaussNoise4 (index) =rand;      %N (0, 1) 
end 
RayleighNoise2=sqrt ( (GaussNoise3)  . ~2  + (GaussNoise4)  . ~2)  ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   Make  Rician  amplitude  of  Channel  1  delayed  pulse  train 
%   and  remove  the  dc  component  of  entire  signal. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SE=0; 

for  index=l :NumberPointsDel,   %1  -->  delay  +  delayed  signal 
if  DelayedPulseTrain (index) ==1   %if  one,  make  Rician 

DelayedPulseTrain (index) =DelayedPulseTrain (index) . *rician (Si 
gAmplitude) / 

SE=  (DelayedPulseTrain (index) )  . "2+SE; 
end       %end  Rician  modification  loop 
end 
SP=SE/NumberPointsDel; 

%%%%%%%   Add  noise  to  Channel  1  Delayed  Pulse  Train  %%%%%%% 
for  index=l :NumberPointsDel,   %1  -->  delay  +  delayed  signal 
if  DelayedPulseTrain (index) ==0  %no  Rayleigh  noise  to  signal 

DelayedPulseTrain (index) =DelayedPulseTrain (index) +RayleighNo 
ise (index) / 

end 
end 

%%%%%  Remove  dc  component  of  Channel  1  Pulse  Train  %%%% 
Chldc=mean (DelayedPulseTrain (1 iNumberPointsDel) ) ; 
DelayedPulseTrain (1 :NumberPointsDel) =DelayedPulseTrain (1 :Num 
berPointsDel) -Chide; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Make  Rician  amplitude  of  Channel  2  delayed  pulse  train 
%  and  remove  dc  component  of  entire  signal. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
SE=0; 

for  index=l :NumberPointsDel2,  %1  -->  delay  +  delayed  signal 
if  DelayedPulseTrain2 (index) ==1       %if  one,  make  Rician 

DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) . *rician ( 
SigAmplitude2) ; 

SE= (DelayedPulseTrain2 (index) ) . A2+SE; 
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end   %end  Rician  modification  loop 
end 
SP=SE/NumberPointsDel2  ; 

%%%%%%   Add  nc  3e  to  Channel  2  Delayed  Pulse  Train  %%%%%%%% 

for  index=l :NumberPointsDel2,  %1  -->  delay  +  delayed  signal 

if  DelayedPulseTrain2 (index) ==0  %no  Rayleigh  noise  to  sig 

DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) +Rayleigh 
Noise2 (index) / 

end 
end 

%%%%%  Remove  dc  component  of  Channel  2  Pulse  Train  %%%% 
Ch2dc=mean (DelayedPulseTrain2 (1 :NumberPointsDel2) ) ; 
DelayedPulseTrain2 (1 :NumberPointsDel2) =DelayedPulseTrain2 (1 : 
NumberPointsDel2) -Ch2dc; 

DelayedPulseTrain=DelayedPulseTrain (1 : 100) ; 
DelayedPulseTrain2=DelayedPulseTrain2 (1 : 100) ; 

RayleighNoise=RayleighNoise (1 : 100)  ; 
RayleighNoise2=RayleighNoise2 (1 : 100) / 
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APPENDIX  C.  TDOA  CONSTANT  THRESHOLD  SIMULATION 

A.    INTRODUCTION 

Four  graphs  of  the  output  of  a  time  domain  based 
correlation  detector  using  a  constant  threshold  are  given. 
Each  figure  has  four  simulations  using  a  SNR  of  ldB.  The  P£a 
is  varied  from  0.01  to  0.00001  from  Figure  26  to  Figure  29. 
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Figure  26.   Time  domain  based  correlation  detector  output 
for  a  SNR  of  ldB  and  a  P£a  of  0.01. 
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Figure  27.   Time  domain  based  correlation  detector  output 
for  a  SNR  of  ldB  and  a  P£a  of  0.001. 
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Figure  28.   Time  domain  based  correlation  detector  output 
for  a  SNR  of  ldB  and  a  Pfa  of  0.0001. 
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Figure  29.   Time  domain  based  correlation  detector  output 
for  a  SNR  of  ldB  and  a  Pfo  of  0.00001. 
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APPENDIX  D.  FFT  OUTPUT  NOISE  TO  SIGNAL  RATIO 

A.   INTRODUCTION 

The  noise  to  signal  ratio  for  thirteen  passes  through  a 
FFT  butterfly  are  calculated  below.  The  calculation  assumes 
the  word  length  of  the  arithmetic  logic  unit  (ALU)  of  the 
computer  is  16  bits  long.  The  length  of  the  input  data  points 
are  12  bits  long.  The  data  points  are  right  justified  when 
placed  into  the  ALU.  For  a  16  bit  ALU  and  a  12  bit  data 
input,  scaling  occurs  at  the  seventh  pass  and  every  other  pass 
after  that. 
Pass  1 

The  noise  power  at  the  output  of  the  butterfly  after  the 
first  pass  NA2(1)  is  a  function  of 

i\TJ(l)  =2WJ(0)  +  2a2s(0)    +  S2aU0)    +  O2r(0) 

where  S2  =  signal  power  input  to  butterfly 

Or  -  truncation  noise  power  (74) 

o2s  =  scaling  noise  power 
o2w  -  sine/ cosine  table  noise  power 
NJiiO)    =  noise  power  input  to  FFT    . 

During  pass  one  there  is  no  scaling,  no  truncation,  and  the 
trig  table  is  not  used.  These  error  terms  are  zero.  Equation 
74    reduces    to 
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ATA2(1)     =   2AT~ 

where  NX  =   ±2'2B  (75) 

Q         3 

=  noise  input  to  butterfly     . 

Pass    2 

The  noise  power  at  the  output  of  the  butterfly  after  the 
second  pass  NA2(2)  is  a  function  of 

Wj(2)  =2iVJ(l)  +  2o|(l)  +  2S202/(D  +  cr£(l)   .      (76) 
During  pass  two  scaling  and  truncation  are  not  performed,  and 
the  trig  table  is  not  used.  These  noise  terms  are  zero. 
Substituting  Equation  75  into  76 


tfj(2)  =  2{2N*\ 


=  222V2 


(77) 


Pass  3 

The  noise  power  at  the  output  of  the  butterfly  after  the 
third  pass  NA2(3)  is  a  function  of 

Nj(3)  =  2WJ(2)  +  2o25(2)  +  22S2al(2)    +  o2r(2)   .      (78) 
Scaling  noise  is  zero  because  it  is  not  performed  during  pass 
three.  Substituting  Equation  77  into  78 
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Nj(3)    =  2{22N2Q)    +  2zS2a2w  +  o2T 

(79) 

=  S3N2   +  22S2a2w  +  o^.     . 


Pass    4 

The  noise  power  at  the  output  of  the  butterfly  after  the 
fourth  pass  NA2(4)  is  a  function  of 

ZVf(4)  =  2WJ(3)  +  2o2(3)  +  23S202,(3)  +  o2T(3)   .      (80) 
Scaling  noise  is  zero  because  it  is  not  performed  during  pass 
four.  Substituting  Equation  79  into  80 

ATJ(4)  =  24N2  +  24S202  +  3o2r  .  (81) 

Pass  5 

The  noise  power  at  the  output  of  the  butterfly  after  the 
fifth  pass  NA2  (5)  is  a  function  of 

WJ(5)  =  2NJ(4)  +  2o2(4)  +  24S202,(4)  +0^(4)   .      (82) 
Scaling  noise  is  zero  because  it  is  not  performed  during  pass 
five.  Substituting  Equation  81  into  82 

N}{5)    =  25N*   +  (25  +  2*)S2o2w   +  la\     .  (83) 

Pass  6 

The  noise  power  at  the  output  of  the  butterfly  after  the 
sixth  pass  NA2(6)  is  a  function  of 
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tfj(6)  =2iVJ(5)  +  2o2(5)  +  2?S2o2/(5)  +  o2r(5)   .      (84) 
Scaling  noise  is  zero  because  it  is  not  performed  during  pass 
six.  Substituting  Equation  83  into  84 

tfj(6)  =  26W2  +  27S2a2w   +  15o2r  .  (85) 

Pass  7 

The  noise  power  at  the  output  of  the  butterfly  after  the 
seventh  pass  NA2(7)  is  a  function  of 

7VJ(7)  =2A£(6)  +  2o2(6)  +  26S2o2/(6)  +0^(6)   .      (86) 
Scaling  is  performed  for  the  first  time  during  pass  seven. 
Substituting  Equation  85  into  86 

WJ(7)  =  27N2   +  (28  +  26)S2o2  +  3102r  +  2o2   .        (87) 

Pass  8 

The  noise  power  at  the  output  of  the  butterfly  after  the 
eighth  pass  NA2(8)  is  a  function  of 

Wj(8)  =  2WJ(7)  +  202s(7)    +  275202,(7)  +0^(7)   .      (88) 
Scaling  is  not  performed  during  this  pass.  Only  the  scaling 
noise  from  the  previous  pass  is  added.  Substituting  Equation 
87  into  88 

Nj(8)  =  2*N2Q   +  (29  +  28)S202  +  63o2r  +  4o2   .        (89) 
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Pass  9 

The  noise  power  at  the  output  of  the  butterfly  after  the 
ninth  pass  NA2(9)  is  a  function  of 

/VJ(9)  =2JVJ(8)  +  2a2s(8)    +  28S2o2„(8)  +  o2r(8)   .      (90) 
Substituting  Equation  89  into  90 

NJ(9)  =  29N2  +  (210  +  29  +  2e)S2o2  +  127o2r  +  10o2s     .    (91) 

Pass  10 

The  noise  power  at  the  output  of  the  butterfly  after  the 
tenth  pass  NA2(10)  is  a  function  of 

wf(10)  =2iVJ(9)  +  2a2s{9)    +  29S2Oz,(9)  +  02r(9)   .     (92) 
Scaling  is  not  performed  this  pass.  Only  the  scaling  noise 
from  the  previous  pass  is  added.  Substituting  Equation  91  into 
92 

tfj(10)  =  2l0N*   +  2l2S2ol   +  25502-  +  20o2g  .        (93) 
Equation  93  describes  the  noise  power  at  the  output  of  a 
butterfly  after  the  tenth  data  pass.  The  noise  to  signal  ratio 
after  the  tenth  pass  through  the  butterfly  is  expressed  as  the 
ratio 
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n!                    21QNt   +  212S2a2w  +  2550?.  +  2Qol  ,QA. 

_J1(10)    =   ? 1 1 £      ,  (94) 

S2  210S2 

Dividing  through  by  the  denominator   (signal  power)   each 
independent  and  uncorrelated  noise  term  can  be  identified 

— 1  (10)  =  — 2  +  Aozw   +  —  —  +  — 

52        S2         4  s2  51  52 

where  —  (10)    =  output  Noise  to  Signal  ratio 

S2 

— -  =  input  Noise  to  Signal  ratio  (95) 

4c2,  =  sine/ cosine  table  noise 

_2 

1    Oj  .  .  .  ,  . 

= truncation  Noise  to  Signal  ratio 

4  s2 

1  o| 
= scaling  Noise  to  Signal  ratio     . 

51  s2 

The  noise  to  signal  ratio  calculation  for  13  passes  (213 
=  8192  points)  through  the  butterfly  algorithm  continues 


Pass  11 

The   noise  power   at   the   output    of  the  butterfly   after  the 
11th  pass   NA2(11)    is   a   function   of 

Nj(ll)    =2NJ(10)    +  2a2s(10)    +  2 10S2O?,(  10)    +  a2r(10)      .       (96) 
Substituting  Equation    93    into    96 

Nj  (11)    =  2X1N2   +    (213   +  210)S2O2r  +  511a2T  +  42a2      .  (97) 
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Pass  12 

The  noise  power  at  the  output  of  the  butterfly  after  the 
12th  pass  NA2(12)  is  a  function  of 

Nl  {12)    =  2NJ(11)  +  2(7*  (11)  +  211S2a2,(ll)  +  a2T(ll)   .   (98) 

Scaling  is  not  performed  this  pass.  Only  scaling  noise  from 
the  previous  pass  is  added.  Substituting  Equation  97  into  98 

NJ(12)    =  212tf2   +    (214   +  212)S2<72f  +   1023<Jt  +   84a|     .  (99) 

Pass    13 

The  noise  power  at  the  output  of  the  butterfly  after  the 
13th  pass  NA2(13)  is  a  function  of 

(100) 

Nj|<13)    =  2Ni(12)    +  2a2s(12)    +  212S2o£(12)     +  a2T(12)       . 

Substituting   Equation    99    into    100 

(101) 

Nj(13)    =  213N2   +    (215   +  213   +  212)52ofr  +  2047a2-  +  170a2,     . 

Equation  101  describes  the  noise  power  at  the  output  of  a 
butterfly  after  the  13th  data  pass.  The  noise  to  signal  ratio 
after   the    13th  pass    is   expressed  as   the   ratio 

Nj  213N2   +    (215  +  213  +  212)52a2r  +  2047a2T  +  170a2;     (102) 

"S1         }  2l3S2 

Dividing  through  by  the  denominator,  each  independent  and 
uncorrelated  noise   term  can  be    identified 
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»t2  »t2  ^-2  ~2 

_ 1  (13)  =  —  +  5.5a2,  +  —  —  + 
S2  S2  4  S2         48.2  S2 

where  —  (13)  =  output  Noise  to  Signal  ratio 

N2 

— 1   =  input  Noise  to  Signal  ratio  (103) 

5 .  5G2,   =  sine/cosine  ta&Ie  noise 

_2 

1  <Jt  .      . 

= truncation  Noise  to  Signal  ratio 

As2 

1        G2s  . 

= scaling  Noise  to  Signal  ratio     . 

48.2  S2 
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